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Preface 


Due  to  the  highly  theoretical  nature  of  the  problem,  this  work  does  not  deal 
with  the  practical  application  of  mixed  Hj/H^  theory  at  all.  I  have  enjoyed 
performing  the  analytical  research,  and  I  trust  the  reader  will  be  able  to  see 
through  the  theory  to  the  possible  future  applications  of  this  particular  discipline. 
I  am  convinced  that  as  the  mixed  theory  becomes  more  mature,  it  will 

provide  large  breakthroughs  in  the  controls  community  and  find  its  way  into  all 
kinds  of  applications.  Since  I  have  been  able  to  contribute  even  a  small  part 
toward  this  end,  this  research  has  been  as  fulfilling  as  it  has  been  trying. 

We,  as  engineers,  spend  all  our  energies  trying  to  find  some  manageable 
model  of  the  phenomena  we  observe  in  nature  with  the  hopes  of  being  able  to 
design  something  useful.  The  more  I  study  the  sciences  (and  in  this  research 
project  in  particular)  the  more  impressed  I  am  with  the  awesome  complexity  of 
the  natural  world  in  which  we  live.  Even  some  of  the  most  elementary  things 
still  remain  a  mystery  to  us.  Since  I  have  been  able  to  pull  back  the  veil  and 
reveal  a  tiny  piece  of  new  information,  I  feel  that  this  thesis  is  more  than  just  a 
technical  work  that  will  sit  on  some  shelf.  It  is  another  chapter  in  the  unfolding 
testimony  of  the  marvelous  God  who  conceived  all  of  this  and  brought  it  into 
being. 

There  have  been  several  people  who  have  been  crucial,  without  whom  this 
thesis  would  not  have  been  possible.  First  is  my  thesis  advisor,  Capt  Ridgely. 

ii 


I  would  like  to  thank  him  for  his  guidance,  patience,  and  technical  expertise. 
This  is  really  his  inspiration,  and  I  am  glad  I  was  able  to  share  in  it.  He  devoted 
more  time  to  his  thesis  students,  just  to  get  us  up  to  speed,  than  any  student 
could  possibly  ask.  Next,  I  extend  my  sincere  thanks  to  my  committee  members, 
Maj  Mracek  and  Dr  Liebst.  Maj  Mracek’s  knowledge  of  the  problem,  numerics, 
and  the  computer  code  is  was  what  got  me  going-and  then  kept  me  going.  His 
personal  interest  in  my  work  is  greatly  appreciated.  Dr  Liebst  sacrificed  many 
hours  and  provided  a  fresh  new  look  at  the  problem  at  a  point  when  it  looked  like 
all  avenues  had  been  exhausted. 

Finally,  I  would  like  to  thank  the  one  who  sacrificed  more  for  this  thesis  than 
anyone,  my  wife  Sherry.  It’s  one  thing  for  a  man  to  have  a  wife,  but  it’s  a 
whole  other  thing  for  him  to  have  a  helpmate  and  friend.  I  appreciate  her 
strength  and  love  more  than  I  can  possibly  say.  I  look  forward  to  being  able  to 
make  up  lost  time  with  her  and  my  special  little  boy,  Ben. 


Scott  R.  Wells 
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Abstract 


0,  rr  r< 


The  problem  of  minimizing  the  2-norm  of  one  transfer  function  subject  to  an 
-a»-norm  bound  on  another  transfer  function  is  e^mined  for  increased  order 
controllers.  In  particular,  the  theoretical  results  of  the  full  order  case  are 


extended  to  the  higher  order  case,  and  SISO y&nd/MIMO  numerical  examples  are 
given  for  increasingly  higher  order  compensators.  Some  of  the  key  proofs  for 
higher  order  compensators  include:  the  global/ minimum  2-norm  is  unachievable 


under  output  feedback  for  certain  levels  of  i  regardless  of  compensator  order; 


the  solution  to  the  mixed  H2/H<^  problem  lies  on  the  boundary  of  the  iso -norm 
constraint  for  this  same  range  of  Y’s;  and  the  suboptimal  mixed  problem 
converges  to  the  optimal  in  the  limit  for  higher  order  controllers.  Also,  it  is 
shown  that  the  optimal  compensator  order  for  the  mixed  H2/H^  problem  is 
greater  than  the  order  of  the  plant  under  certain  conditions,  and  a  conjecture 


about  the  optimal  order  for  the  mixed  problem  is  made. 
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INVESTIGATION  OF  THE  EFFECTS  OF 
INCREASED  ORDER  COMPENSATORS  IN 
MIXED  H^H*,  OPTIMIZATION 


I.  Introduction 


1.1  Background 

Optimal  control  theory  provides  powerful  tools  for  designing  feedback 
controllers,  particularly  when  the  dynamic  system  being  controlled  is  very 
complex.  Classical  methods  such  as  root  locus  become  too  cumbersome  or 
completely  unusable  for  high  order,  multiple-input  multiple-output  systems,  while 
optimization  techniques  can  offer  a  mathematical  structure  that  readily  handles 
these  complicated  systems.  In  an  optimal  control  synthesis  problem,  a  controller 
that  will  minimize  some  prescribed  cost  objective  (sometimes  with  additional 
constraints)  is  sought.  Extreme  care  must  be  taken  when  defining  this  cost 
function  because  any  controller,  including  one  that  results  in  unacceptable  time 
responses  (or  even  an  unstable  closed-loop  system),  can  be  shown  to  be  optimal 
to  some  cost  function.  Obviously,  the  key  in  optimal  control  is  defining  the 
"right"  cost  objective.  Two  performance  measures  which  are  currently  receiving 
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a  great  deal  of  attention  in  the  controls  community  are  the  H2  and  H^,  norms. 
These  norms  are  defined  by 


li  G(ju)  II 2  s  [  °°  tr  [G  *  (j“)  G(jw)]  do 

Z7T  J  -<*> 

|G<j«)l„  -  uSgPR  5  (G(i«)] 

(  a  =  max  singular  value ) 


[Dai90] 


Both  of  these  performance  measures  are  well  motivated  and  have  significant 
merit.  By  themselves,  however,  each  produce  controllers  that  have  potentially 
undesirable  characteristics.  This  has  led  to  the  development  of  the  mixed 
optimization  problem. 

While  it  is  beyond  the  scope  of  this  work  to  review  the  vast  amount  of 
literature  on  H2  and  theory,  some  of  the  key  ideas  need  to  be  summarized 
in  order  to  properly  motivate  the  mixed  problem.  Consider  the  general  feedback 
control  system  shown  in  Figure  1-1.  P  is  the  generalized  linear,  time-invariant 
plant  transfer  function  and  contains  all  the  weighting  functions  required  to  obtain 
this  general  form.  The  exogenous  and  controlled  input  vectors  are  w  and  u, 
respectively.  z  and  y  are  the  controlled  and  measured  output  vectors, 
respectively.  Denote  the  closed-loop  transfer  function  from  w  to  z  as  Tzw. 
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z 


y 


Figure  1-1.  H2  Feedback  Control  System  Block  Diagram 


If  w  is  assumed  to  be  white  Gaussian  noise  and  certain  assumptions  are  made  on 
the  plant,  the  H2  optimization  problem  can  be  formulated  and  solved.  The  cost 
function  here  is  the  2-norm  of  the  closed-loop  transfer  function  Tzw.  The  only 
compensators  that  are  considered  "admissible"  in  the  optimization  are  those  that 
are  real-rational,  proper,  and  internally  stabilizing.  The  problem,  stated  more 
formally,  is:  infimize  the  2-norm  of  Tzw  over  the  set  of  all  admissible 
compensators,  or  find  the  K  that  achieves 

inf  ||  Tzw  1 2  ■  a0 
K  adm 


Note  from  the  definition  of  the  2-norm  that,  given  a  unit  intensity  white  noise 
input,  (by  Parseval’s  theorem)  the  square  of  ||Tzw||2  is  equal  to  the  energy  of  the 
output  signal.  The  H2  optimal  controller  is  the  one  that  results  in  the  minimum 
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energy  of  the  controlled  output  z  due  to  the  input  w.  This  is  very  desirable  from 
a  performance  point  of  view  for  white  noise  input  applications. 

It  can  be  shown  [Rid91a,206-216]  that  H2  optimization,  under  the  assumption 
of  output  feedback,  is  equivalent  to  a  corresponding  Linear  Quadratic  Gaussian 
(LQG)  problem.  That  is,  the  compensator  that  minimizes  the  LQG  cost  function 

J  =  jo°cCrrQ;t  +  utRu  +  *tNk)  dt 

is  an  H2  optimal  compensator.  Both  problems  involve  the  solution  of  two  Riccati 
equations  and  produce  a  unique  controller  whose  order  is  equal  to  the  order  of 
the  plant  (full  order).  It  is  well  known  that  while  Linear  Quadratic  Regulators 
(LQR)  and  Estimators  (LQE)  generate  systems  with  guaranteed  minimum  gain  and 
phase  margins,  LQG  compensators  do  not  exhibit  this  same  quality.  In  fact,  it 
is  possible  for  an  LQG  system  to  have  gain  and  phase  margins  that  are  arbitrarily 
small  [Doy78, 756-757].  So,  while  the  H2  optimal  controller  has  desirable 
performance  characteristics,  it  does  not  provide  any  guarantees  on  system 
robustness. 

It  is  possible  in  the  LQG  problem  to  recover  to  some  degree  the  guaranteed 
margins  of  the  LQR/LQE  systems  by  using  a  technique  called  Loop  Transfer 
Recovery  (LTR).  In  LTR,  the  stability  margins  are  recovered  by  trading  off 
regulator  or  estimator  performance  (depending  on  where  the  model  uncertainties 
are  entering  the  system).  In  LQG/LTR  the  designer  must  decide  where  to  break 
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the  loop  (at  the  input  of  the  plant  or  at  the  output)  when  defining  measures  of 
performance  and  stability  margins.  Unfortunately,  it  is  not  possible  to  define 
performance  at  one  location  and  stability  margins  at  another  [RB86,9. 1-9.8].  So, 
while  LQG/LTR  is  an  effective  design  method  for  certain  cases,  it  cannot  handle 
the  most  general  case. 

Consider  now  the  problem  of  optimization.  The  generalized  block 
diagram  is  shown  in  Figure  1-2. 


e 


y 


Figure  1-2.  Feedback  Control  System  Block  Diagram 


The  exogenous  input  is  defined  to  be  the  vector  d,  which  is  assumed  to  be  a 
bounded  energy  input.  The  controlled  output  vector  is  denoted  by  e.  The 
problem  set-up  is  exactly  the  same  as  before,  except  now  it  is  the  oo-norm  of  the 
closed-loop  transfer  function  Ted  that  is  being  minimized;  that  is  find  a  K  which 
achieves 

inf  II  Ted  ||  a,  =  70 
K  adm 
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The  true  solution  to  this  optimization  problem  is  usually  avoided,  and  there  are 
several  reasons  why.  First,  the  design  algorithms  are  necessarily  iterative, 
and  achieving  the  exact  optimal  value  is  very  difficult.  Also,  optimal 
controllers  tend  to  have  some  undesirable  characteristics.  They  are  typically 
infinite  bandwidth  compensators  which  produce  all-pass  closed-loop  maximum 
singular  value  plots  [ZDGB90,  2502-2503].  Thus,  the  true  optimal  solution 
is  not  only  difficult  to  calculate,  but  it  may  be  impractical  for  a  real  system.  The 
more  practical  solution  is  the  H^,  swAoptimal  controller,  which  is  the 
compensator  that  insures 

II Ted  II  oo  <  7,  where  y  >  y0 

Both  the  suboptimal  and  optimal  compensators  are,  in  general,  non-unique. 
There  are  an  infinite  number  of  controllers  that  achieve  a  given  level  of 
performance. 

One  of  the  key  advantages  of  the  oo-norm  is  its  direct  link  to  system 
robustness.  This  is  possible  because  the  oo-norm  has  the  submultiplicative 
property.  That  is,  for  some  F,G  E  RL*,,  where  RL^  is  the  Banach  space  of 
all  real-rational  proper  stable  transfer  matrices, 

|| FOIL  £  II F ||  0.  || G ||  g.  [Zam66] 
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Note  that  the  2-norm  does  not  have  this  property.  To  see  the  application  of  this 
property,  consider  the  uncertainty  block  diagram  shown  in  Figure  1-3.  Ted  is  the 
nominal  closed-loop  transfer  function,  and  A  is  a  perturbation  function  that 
characterizes  the  system’s  uncertainty. 


d 


A  k 


j  a 

1 

1 

1  I 

I  Ted  | 

1 _ 1 

e 


Figure  1-3.  Uncertainty  Block  Diagram 

The  Small  Gain  Theorem  gives  a  relationship  that  describes  how  large  A  can  be 
before  the  nominal  system  becomes  unstable. 

Small  Gain  Theorem:  Assume  Ted(s),  A(s)  G  RH*. 

If  ||  Ted(s)  A(s)  |L  <  1 

then  the  closed-loop  system  (with  A)  is  stable  [Zam66] 

Then,  by  the  submultiplicative  property  of  the  oo-norm, 

II  Ted(s)  A(s)|L  £  l|Te,(s)|L  |A(s)|L  <  1 
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Thus,  to  ensure  system  stability, 

II  A(s)  ||  «  <  _L_ 

II  Ted  II  oc 

The  smaller  ||Tedj|  «*,  is,  the  more  uncertainty  the  system  can  handle.  Due  to  the 
submultiplicative  property  of  the  oo-norm,  optimization  seems  like  an 
obvious  choice  for  designing  a  system  for  stability  robustness. 

As  can  be  seen,  two  reasonable  choices  for  the  cost  objective  in  the 
optimization  problem  are  the  2-norm  and  co-norm  of  the  closed-loop  transfer 
functions.  The  2-norm  is  desirable  because  it  results  in  optimal  (in  the  energy 
sense)  performance  in  the  face  of  white  noises.  The  oo-norm  is  desirable  because 
it  guarantees  a  level  of  robustness  to  plant  uncertainties.  What  would  be  most 
desirable  would  be  a  methodology  that  incorporates  both.  This  is  what  has  been 
termed  the  mixed  optimization  problem.  Consider  now  the  general  form 

of  the  feedback  control  system  for  the  mixed  problem  as  shown  in  Figure  1-4. 


Figure  1-4.  Mixed  ftyH*  System  Block  Diagram 
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The  exogenous  input  d  is  assumed  to  be  a  deterministic  signal  of  unknown  but 
bounded  energy.  The  input  w  is  assumed  to  be  a  zero-mean  white  Gaussian  noise 
of  unit  intensity.  The  signals  e  and  z  are  controlled  outputs,  and  may  be  equal, 
dependent,  or  independent.  The  signals  u  and  y  are  the  control  inputs  and 
measured  outputs,  respectively.  All  weighting  functions  required  to  obtain  this 
general  form  are  included  in  the  generalized  plant  transfer  function  P.  Denote 
the  closed-loop  transfer  functions  from  w  to  z  and  from  d  to  e  as  Tzw  and  Ted, 
respectively.  Stated  formally,  the  general  mixed  H2/Hao  synthesis  problem  is  to 
find  an  admissible  controller  K(s)  that  achieves 

inf  ||  Tzw  ||  2,  subject  to  the  constraint  ||  Ted  ||  »  <  7 
K  adm 

As  will  be  shown  later,  in  the  region  of  interest  ||TZW||2  and  ||  Ted  ||  ^  are 
competing  objectives.  This  makes  sense  physically  since  it  is  expected  that  some 
performance  should  have  to  be  sacrificed  in  order  to  gain  robustness  (and  vice 
versa).  The  mixed  H2/Hoc  problem  provides  the  needed  structure  for  the 
designer  to  explicitly  observe  and  influence  this  trade-off.  An  interesting  note 
is  that  if  the  problem  is  set  up  such  that  d  =  w  and  e  =  z,  the  closed-loop 
transfer  functions  Ted  and  Tzw  are  equal.  The  problem  is  then  essentially 
equivalent  to  LQG/LTR,  and  the  recovery  of  the  stability  margins  is  seen  directly 
in  the  trade-off  between  the  2-norm  and  00 -norm.  The  nice  thing  about  the 
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Hj/Hg,,  approach  is  that  the  input  and  output  signals  can  be  defined  in  any  way 
the  designer  chooses,  so  the  limitations  of  LQG/LTR  are  not  present  and  the 
problem  remains  completely  general. 

The  value  of  designing  optimal  controllers  with  mixed  H2  and  H^, 
performance  objectives  has  been  recognized  for  some  time,  but  the  solution  to  the 
most  general  problem  had  not  been  demonstrated  until  recently  by  Ridgely  in 
[Rid91a].  Due  to  the  complexity  of  the  problem,  his  solution  requires  numerical 
techniques.  In  fact,  it  is  generally  believed  that  the  problem  may  not  have  an 
explicit  analytical  solution. 

Ridgely  was  not  the  first  to  examine  the  mixed  problem;  several 

earlier  papers  appeared  starting  in  1989.  The  early  works  by  Bernstein  & 
Haddad  ([BH89]);  Zhou,  Doyle,  Glover  &  Bodenheimer  ([ZDGB90]);  Yeh,  Banda 
&  Chang  ([YBC90]);  Mustafa  &  Glover  ([MG90]);  and  Khargonekar  &  Rotea 
([KR91])  laid  the  foundations  for  the  general  mixed  problem,  but  made  varying 
assumptions  that  specialized  the  more  general  problem.  In  addition,  they  did  not 
minimize  the  actual  2-norm  but  rather  an  upper  bound  to  it.  It  was  not  until 
Rotea  &  Khargonekar’s  work  ([RK9 1  ])  that  a  special  case  of  the  true 
nonconservative  problem  was  solved.  They  allowed  two  sets  of  inputs  and 
outputs  and  did  not  use  an  upper  bound  to  the  2-norm;  however,  they  did  restrict 
their  solution  to  the  case  of  full  state  availability  for  feedback.  Finally, 
Ridgely’s  work  offered  the  first  solutions  to  the  general  nonconservative  problem. 
He  used  a  Lagrange  multiplier  technique  and  derived  the  set  of  necessary 
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conditions  for  the  output  feedback  case.  He  also  proved  some  key  properties  for 
full  order  compensators  and  performed  an  exhaustive  analysis  of  a  SISO  and 
MIMO  example,  providing  valuable  insight  into  the  nature  of  the  problem. 

While  the  bulk  of  Ridgely’s  work  was  done  with  full  order  controllers,  he  did 
take  a  cursory  look  at  higher  order  compensators.  He  found,  by  example,  a 
higher  order  controller  that  satisfied  the  necessary  conditions  and  produced  a 
2-norm  lower  than  the  corresponding  full  order  compensator.  It  is  obvious  from 
this  example  and  other  analysis  that  the  optimal  order  of  the  mixed  compensator 
(that  is,  the  lowest  order  compensator  with  a  minimal  realization  that  achieves 
the  true  optimal  mixed  solution  2-norm)  is  not  the  order  of  the  plant.  However, 
the  ultimate  issue  of  optimal  order  was  not  determined. 

1.2  Research  Objectives 

The  work  contained  here  is  an  extension  of  Ridgely’s  dissertation.  His 
solution  methodology  was  used  and  was  applied  specifically  to  compensators 
whos!  order  is  higher  than  the  plant.  The  main  objective  of  the  research  was  to 
use  the  example  systems  defined  in  Ridgely’s  dissertation  and  numerically  solve 
the  mixed  H2/HO0  for  increasingly  greater  order  controllers.  An  investigation  of 
what  these  results  say  about  the  nature  of  the  problem  and  what  the  optimal  order 
might  be  was  then  performed.  Obviously,  these  numerical  solutions  cannot 
directly  prove  optimal  order  in  general,  but  the  investigation  could  shed  new  light 
on  this  difficult  subject.  In  addition  to  running  numerical  examples,  another 
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objective  was  to  extend  some  of  the  analytical  proofs  for  the  full  order  case  to 
the  higher  order  case.  Finally,  while  it  was  not  expected  to  be  able  to  formally 
prove  the  optimal  order  of  the  mixed  problem,  this  question  motivated  all 

the  research,  and  the  proof  (or  at  least  a  strong  conjecture)  was  always  an 
underlying  objective. 

1.3  Thesis  Outline 

This  thesis  is  comprised  of  six  chapters.  Chapter  I  provides  the  background 
for  this  work  and  outlines  the  research  accomplished.  The  motivation  for  doing 
mixed  optimization  is  presented.  Also,  a  brief  history  of  the  problem 

development  is  given. 

Chapter  II  gives  motivations  for  investigating  higher  order  compensators  and 
takes  an  introductory  look  at  the  question  of  optimal  order.  Some  of  the 
problems  with  determining  optimal  order  are  discussed.  Then,  the  optimal  orders 
of  related  optimization  problems  are  reviewed.  Specifically,  the  optimal  orders 
of  compensators  in  the  pure  H2,  Hw,  and  minimum  entropy  problems  are  shown. 

Chapter  III  then  begins  the  formal  development  of  the  mixed  problem.  After 
defining  the  problem  statement,  first-order  necessary  conditions  for  the  general 
and  suboptimal  mixed  problems  are  derived.  Then  a  brief  review  of  the  full 
order  analytical  results  is  given. 

Chapter  IV  contains  all  the  analytical  proofs  that  extend  the  full  order  results 
to  higher  order  compensators.  The  true  optimal  mixed  problem  is  the  main  focus 
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of  theory  presented  here;  however,  the  suboptimal  mixed  problem  is  also  briefly 
discussed. 

Chapter  V  is  the  section  that  contains  all  the  numerical  results  of  the 
research.  First,  the  algorithm  used  to  numerically  solve  the  set  of  necessary 
conditions  is  discussed.  Then,  a  SISO  and  MIMO  example  are  given  along  with 
discussions  of  their  results. 

Finally,  Chapter  VI  concludes  the  work  with  a  formal  conjecture  of  optimal 
order  supported  by  both  analytical  and  numerical  evidence.  Recommendations 
for  future  research  are  offered  and  a  summary  is  given.  The  FORTRAN  source 
code  of  the  numerical  algorithm  discussed  in  Chapter  V  is  included  as  an 
appendix. 
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II.  Higher  Order  Compensators  and 
Optimal  Order 


2.1  General  Discussion 

Why  even  consider  higher  order  compensators?  From  a  practical  point  of 
view,  it  is  unlikely  that  a  higher  order  controller  would  ever  be  implemented  on 
a  real  system,  especially  if  the  plant  has  any  significant  size  at  all.  For  example, 
consider  a  typical  aerospace  application,  aircraft  pitch  attitude  control.  The 
nominal  plant  from  the  linearized  perturbation  equations  is  fourth  order.  After 
adding  in  all  the  weights  for  the  input  and  output  signals,  the  general  system 
plant  could  easily  double.  Now,  the  full  order  controller  must  be  8th  order. 
This  is  already  becoming  unreasonable  from  a  real-life  application  point  of  view. 
Now,  say  that  the  optimal  order  of  a  controller  for  some  optimization  problem 
is  three  times  the  order  of  the  plant.  This  would  be  a  24th  order  controller.  If 
it  is  even  possible  to  build  such  a  compensator,  it  would  undoubtedly  be 
prohibitively  expensive.  It  seems  the  more  useful  area  of  interest  for  practical 
applications  would  be  reduced  order  controllers.  So  again,  why  consider  higher 
order  compensators? 

There  are  at  least  two  very  good  reasons  to  examine  higher  order 
compensators.  The  first  reason  is  more  philosophical  than  practical.  What  is  the 
optimal  order  of  the  compensator  for  the  mixed  H2/H0B  control  problem?  It  has 
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already  been  shown  that  it  is  something  greater  than  the  order  of  the  plant,  in 
general.  Therefore,  if  the  question  of  optimal  order  is  to  be  addressed  (which 
it  must  be--purely  for  the  sake  of  fully  developing  theory),  the  higher 

order  case  must  be  examined. 

The  second  reason  is  more  practical.  The  design  process  is  really  an  art  of 
compromise.  The  engineer  strives  to  design  a  system  that  will  meet  some  set  of 
specifications  without  over-exceeding  the  specifications.  Why  design  a  system 
that  can  handle  ten-fold  variations  in  the  system  parameters  when  three-fold 
variations  are  all  that  are  actually  expected?  The  main  reason,  besides  cost,  for 
not  overdesigning  a  system  is  that  specifications  are  almost  always  competing. 
For  example,  performance  and  robustness  levels  are  two  specifications  that  a 
control  system  would  have  to  satisfy.  Unfortunately,  in  order  to  get  one,  there 
usually  will  have  to  be  sacrifices  made  in  the  other.  In  order  for  an  engineer  to 
properly  make  these  trade-offs,  he  needs  to  know  the  limits  for  the  problem,  such 
as  the  maximum  level  of  performance  and  the  maximum  level  of  robustness  that 
can  be  achieved  by  any  controller.  With  these  parameters  in  hand,  the  engineer 
can  then  back  off  on  both  ends  until  a  compromise  is  found  that  will  satisfy  both 
specifications  (or  determine  that  the  problem  has  been  overspecified).  Thus, 
even  if  the  optimal  order  of  the  mixed  problem  turns  out  to  be  infinite  (which 
could  never  be  implemented  on  a  real  system),  this  infinite  order  compensator 
represents  the  limit  of  achievable  performance  for  a  given  level  of  robustness  and 
is  important  to  know.  Even  more  important  is  the  trade-off  of  performance 
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versus  compensator  order  from  full  order  to  the  optimal  order.  This  would  give 
the  designer  the  ability  to  explicitly  see  what  advantages  can  be  gained  and  how 
much  it  will  cost  to  attain  them.  Therefore,  since  the  optimal  order  appears  to 
be  greater  than  the  order  of  the  plant,  increased  order  compensators  must  be 
investigated,  even  for  practical  design  purposes. 

Determining  the  optimal  order  of  a  compensator  that  satisfies  some 
optimization  criterion  is  not  an  easy  or  straightforward  task.  If  the  problem  can 
be  posed  such  that  compensator  order  does  not  have  to  specified  at  the  beginning 
of  the  problem,  and  the  solution  ends  up  defining  the  compensator  order,  the 
optimal  order  can  be  determined.  However,  if  a  method  like  Lagrange 
multipliers  is  used,  the  order  of  the  compensator  must  be  chosen  during  the  set¬ 
up  of  the  problem.  This  has  the  advantage  of  allowing  the  engineer  to  specify 
the  order  of  the  controller  in  advance  (which  is  particularly  helpful  in  the 
reduced  order  problem),  but  it  completely  prevents  the  determination  of  optimal 
order.  If  one  tries  to  add  compensator  order  as  a  constraint,  a  number  of 
problems  immediately  arise.  First  of  all,  how  does  one  express  order  as  a 
mathematical  relationship  that  can  be  used  in  an  optimization  algorithm?  Also, 
the  sizes  of  the  matrices  in  the  Lagrangian  and  thus  the  necessary  conditions 
change  with  changes  in  order.  This  gives  rise  to  difficulties  in  deriving 
completely  general  proofs.  For  example,  proofs  for  full  order  compensators  that 
require  some  matrix  to  be  square  and  full  rank  may  not  be  valid  for  higher  order 
compensators  because  these  same  matrices  may  not  even  be  square  for  the  higher 
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order  case.  Finally,  compensator  order  is  not  a  convex  constraint.  The  mixed 
H^/Hqo  problem  has  a  convex  objective  function  subject  to  a  convex  constraint 
[Rid91a,27].  However,  if  compensator  order  is  added  as  a  constraint,  the  whole 
nature  of  the  problem  would  be  changed  due  to  this  nonconvex  constraint. 
Therefore,  determining  the  optimal  order  must  be  accomplished  by  discovering 
some  clever  parameterization  of  solutions,  completely  recasting  the  problem  in 
some  different  space,  or  some  form  of  variational  approach. 

In  order  to  provide  insight  into  the  mixed  problem,  first  consider  the  optimal 
order  for  the  related  problems  of  H2,  H^,,  and  entropy  minimization.  The 
discussion  that  follows  is  intended  to  be  more  heuristic  than  rigorous  —  the 
purpose  is  simply  to  provide  the  necessary  background  of  established  theory  and 
to  lay  the  foundation  for  later  discussions  of  the  mixed  problem. 

2.2  H2  Optimization 

This  section  is  divided  into  two  parts.  The  first  part  assumes  availability  of  the 
full  state  vector  for  feedback;  the  other  part  assumes  output  feedback.  The 
nature  of  the  solutions  turns  out  to  be  surprisingly  different.  Both  cases  begin 
with  the  same  system  definitions.  Consider  the  H2  optimization  block  diagram 
given  in  Figure  2-1. 
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Figure  2-1.  H2  Feedback  Control  System  Block  Diagram 

Aga 

in,  the  exogenous  input  w  is  assumed  to  be  a  zero-mean  white  Gaussian  noise  of 
unit  intensity.  The  vector  z  is  the  controlled  output.  The  signals  u  and  y  are  the 
control  inputs  and  measured  outputs,  respectively.  The  state  space  equations  for 
the  generalized  transfer  matrix  P  are  given  by 


dx/dt  =  Ax  +  Bwvv  +Buw 

(2.1a) 

z  =  CjX  +  Dzwvv  -1-  Dzu« 

(2.1b) 

y  =  CyX  +  Dywvv  +  Dyu« 

(2.1c) 

The  transfer  function  matrix  for  the  open-loop  plant  P  can  be  partitioned  as 
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The  closed-loop  transfer  function  from  w  to  z  is  given  by  the  LFT 


T 


zw 


=  Pzw  +  P„  K  [I-PyuK]-'  Pyw 


2.2.1  State  Feedback  Case.  The  solution  is  found  through  a 
parameterization  of  the  set  of  all  admissible  unconstrained  H2-optimal  controllers. 
That  is,  determine  the  set  of  all  K(s)  such  that 

|T2W(K)||2  -  a0 


Assume  the  following  conditions  on  the  plant  P : 


(i)  Cy  =  I,  Dyw  =  Dyu  =  0 
(ii)  (A,BU)  stabilizable 


(iii)  Dzw  =  0 

(iv)  D^UDZU  full  rank 


A-jwI  Bu 

cz  d2U 


(state  feedback) 


full  column  rank  for  all  w  6  R 
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Let  K(s)  denote  an  admissible  controller.  K(s)  can  be  parameterized  by  J  and  the 
freedom  parameter  Q  as  shown  in  Figure  2-2. 


Figure  2-2.  Parameterization  of  K(s) 


From  Theorem  1  of  [RK91],  the  complete  parameterization  of  optimal 
controllers,  K(s),  that  minimize  ||Tzw||2  are  given  by  the  linear  fractional 
transformation  of  J  and  Q  with 


«”  “1 


AF 

0 

Bu 

0 

F 

I 

-I 

I 

0 

and  Q  6  S,  where  S  a  {Q  6  RHa:  Q  =  Wn,(sI-AF),  W  €  RH2} 


F  ^  -(D^D^)-1  (DTZUCZ  +  B l  X) 
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and  where  X  is  the  unique  stabilizing  solution  to  the  ARE: 


ATX  +  XA  -  (DT2UCZ  +  B l  X^D^D^-^D^C,  +  bJ  X)  +  Cz  =  0 

and 

IIj  s  I  -  BUBU+  (where  Bu+  is  the  Moore-Penrose 

pseudoinverse  of  Bu) 

AF  »  A  +  BUF 

Notice  that  when  Im(Bu)  =  Rn  (or  in  other  words,  Bu£Rnxm,  m>n),  this 
parameterization  reduces  to  a  unique,  state-feedback  (static  gain)  controller, 
K0  =  F.  Also,  an  important  result  from  this  parameterization  is  that  when  Im(Bu) 
is  a  proper  subset  of  Rn  (m  <  n),  there  is  a  family  of  (dynamic)  controllers 
parameterized  by  W.  Thus,  the  optimal  order  under  full  state  availability  is  zero. 
Higher  order  compensators  are  part  of  the  set  of  solutions,  but  they  achieve  the 
exact  same  minimum  2-norm  as  the  static  compensator.  If  the  unconstrained 
problem  is  the  only  item  of  interest,  there  would  be  no  reason  for  a  designer  to 
choose  a  higher-order  dynamic  controller  over  the  static  one.  However,  there  is 
extra  freedom  provided  by  this  family  of  H2-optimal  compensators  which  can  be 
exploited  to  satisfy  some  additional  constraints.  This  freedom  leads  to  interesting 
results  in  the  mixed  problem. 

Consider  briefly  the  mixed  fy/H*  problem  (still  under  full  state 
availability).  Rotea  &  Khargonekar  define  two  problems  [RK91, 307-308]. 
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Problem  A:  The  true  mixed  H2/Hoe  problem  —  minimize  the  H2-norm,  subject 
to  an  Ha,  constraint.  That  is,  find 

a*  =  {  inf  ||Tzw||2,  subj  to  ||  Ted  ||  <  7} 

K  adm 

Problem  B;  An  H2  super-optimization  problem  —  determine  the  unconstrained 
H2-optimal  controller  that  also  achieves  an  (added)  H^,  bound.  That  is,  find  all 
K(s)  such  that 

«o  s  {  inf  II  Tzw  II  2>; 

K  adm 

and  the  constraint  ||  Ted  |f  ^  <  7  is  also  trivially  satisfied. 

Note  that  a  solution  to  Problem  B  is  also  a  solution  to  Problem  A  (but  the  reverse 
is  not  necessarily  true). 

Rotea  &  Khargonekar  then  provide  the  necessary  and  sufficient  conditions  for 
the  existence  of  solutions  to  Problem  B  (and  thus  sufficient  conditions  for 
existence  of  solutions  to  Problem  A).  Also,  under  certain  conditions,  existence 
of  a  solution  to  Problem  B  is  necessary  and  sufficient  for  existence  of  a  solution 
to  Problem  A.  If  a  solution  to  Problem  B  exists,  they  show  that  even  though  the 
full  state  is  available  for  feedback,  the  solution  may  necessarily  be  dynamic. 
This  leads  them  to  the  conjecture  that  the  optimal  order  under  output  feedback 
is  greater  than  the  order  of  the  plant. 
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Another  important  thing  to  notice  in  Problem  B  is  that  the  global  minimum 
of  || Tzw |t 2  (a0)  is. achieved  at  each  level  of  j|Ted||  (if  the  solution  exists).  As 
will  be  shown,  the  character  of  the  solution  under  output  feedback  is  considerably 
different. 


2,2.2  Output  Feedback  Case.  First,  parameterize  the  set  of  all  admissible 
unconstrained  H2-.sw£optimal  controllers.  That  is,  determine  the  set  of  all  K(s) 
such  that 

II  Tzw(K)  il  2  ^  ^  ao 

Assume  the  following  conditions  on  the  plant  P: 

(i)  DZw  =  0 

(ii)  Dyu  =  0 

(iii)  (A,BU)  stabilizable,  (Cy,A)  detectable 

(iv)  DTzuDzu,  DywDyw  full  rank 

Without  loss  of  generality,  strengthen  this  so  that 
DTzuDzu  =  I,  DywDjw  =  I 


A  — j<*>  I  Bu 

C2  ^ZU 


full  column  rank  for  all  u  G  R 
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A-jwI  Bw 

^•y  ^yw 


full  row  rank  for  all  u  G  R 


Note  that  the  case  of  full  state  availability  is  not  a  trivial  special  case  of  output 
feedback  for  this  development  since  Dyw  =  0  violates  assumption  (iv)  and  leads 
to  a  singular  control  problem. 

From  [DGKF89],  the  complete  parameterization  of  all  suboptimal  admissible 
controllers,  K(s),  that  achieve  J1  Tzw  (|  2  ^  a,  a  >  aQ,  is  given  by  the  lower  LFT 
of  J  and  Q  shown  in  Figure  2-2,  where 


v~ 


A  J 

Kf 

1 

Kfi 

j  *  ;  K 

!  -  Iv  c 

| 

0 

I 

j_  Kc, 

I 

1 

0 

Aj  —  A  -  Kj-Cy  -  BUKC 

Kc  =  bJ  x2  +  b^c„  Kci  =  -cy 


Kf  =  Y2c^  +  BwdJw, 


Kn  =  Bu 


and  where  X2  and  Y2  are  the  real,  unique,  symmetric  positive  semidefinite 
solutions  to  the  ARE’s: 
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(A-BUDT2UC2)TX2  +  X2(A-B„DT2UC2)  -  x2bubJ  x2 

+  [(I-D,UDT211)CJT[(I-DIUDT211)CJ  =  0 

(A-BwDyTwCy)Y2  +  Y2(A-BwDyTwCy)T  -  Y,Cy  CyY, 

+  [Bw(I-DywDyw)][Bw(I-DywDyw)]T  =  0 
and 

0  6  RH2,  II  Q II I  <  a2  -  ac2 

If  the  H2  optimal  compensator  is  desired  (i.e.  a  =  aG),  then  Q  s  0,  and  the 
solution  is  unique  and  has  order  equal  to  the  order  of  the  plant.  Note  that  this 
is  in  contrast  to  the  state  feedback  case,  where  a  family  of  optimal  compensators 
exists.  If  an  H2  suboptimal  compensator  is  desired  (i.e.  a  >  aQ)  in  output 
feedback,  then  Q(s)  is  not  equal  to  zero,  and  a  family  of  suboptimal  compensators 
exists.  In  summary,  under  the  assumption  of  output  feedback,  the  optimal 
compensator  is  unique  and  the  optimal  order  is  the  order  of  the  plant. 

Now  consider  Rotea  &  Khargonekar’s  Problem  B  for  the  output  feedback 
case.  By  definition  of  the  problem,  a  =  aQ.  Therefore,  QsO  and  K(s)sK2opt 
(a  unique  solution).  Define  the  w-norm  of  Ted  when  the  H2  optimal  compensator 
is  used  as  |  Ted(K2opt)  ||  w  =  72.  If  a  7  is  chosen  such  that  y  <  72*  no  solution 
to  Problem  B  exists  at  all.  If  y  is  chosen  such  that  7  >  y2,  then  K(s)  is  K2opt 
(the  unconstrained  H2-optimal  solution).  While  the  family  of  dynamic  H2  optimal 
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compensators  with  full  state  availability  provides  degrees  of  freedom  that  enable 
7  to  be  reduced  in  Froblem  B,  these  results  are  not  attainable  in  the  output 
feedback  case  because  the  family  of  H2  optimal  compensators  includes  only  one 
unique  controller.  If  Problem  B  has  a  solution,  it  is  unique  with  7*=72. 
Otherwise,  there  is  no  solution. 

2.3  H,0  Optimization 

Consider  the  HB  optimization  block  diagram  given  in  Figure  2-3. 


Figure  2-3.  Feedback  Control  System  Block  Diagram 


Again,  the  exogenous  input  d  is  assumed  to  be  a  deterministic  signal  of 
unknown  but  bounded  energy.  The  vector  e  is  the  controlled  output.  The  signals 
u  and  y  are  the  control  inputs  and  measured  outputs,  respectively.  The  state 
space  equations  for  the  generalized  transfer  matrix  P  are  given  by 
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dr/dt  =  Ax  +  Bdd  +Buw 

(2.2a) 

e  =  CgX  +  D  sdd  +  Deuw 

(2.2b) 

y  =  Cy*  +  Dydt/  +  Dyuu 

(2.2c) 

The  transfer  function  matrix  for  the  open-loop  plant  P  can  be  partitioned  as 


P  = 


The  closed-loop  transfer  function  from  d  to  e  is  given  by  the  LFT 

Ted  =  Ped  +  Peu  K  [I'P^K]'1  Pyd 

Assume  the  following  conditions  on  the  plant  P : 

(i)  Ded  =  0 

(ii)  Dyu  =  0 

(iii)  (A,BU)  stabilizable,  (Cy,A)  detectable 

(iv)  D^uDeu,  DydDyd  full  rank 

Without  loss  of  generality,  strengthen  this  so  that 

Deu^eu  =  I.  DydDyd  =  I 
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A  —  jo>  I  Bu 

®eu 


full  column  rank  for  all  w  €  R 


(vi) 


A-jul 


full  row  rank  for  all  u  €  R 


A  parameterization  of  all  5n6optimal  compensators  will  now  be  given.  As 
discussed  earlier,  the  H^,  optimal  solution  is  not  desired  and  will  not  be  directly 
addressed  here.  The  problem  is  to  find  the  family  of  admissible  compensators 
K(s)  such  that 

II  Ted  |  *  <  7>  where  7  >  70 


From  [DGKF89],  the  complete  parameterization  of  suboptimal  H,*  controllers  is 
given  by  an  LFT  of  J  and  Q  as  shown  in  the  block  diagram  in  Figure  2-2  with 


Aj 

Kf 

Kfi 

J  - 

-Kc 

0 

I 

Kcl 

I 

0 
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where 


A,  =  A  -  K,Cy  -  BUKC  +  -r'2Y„cI  (Ce  -  DeuKc) 

Kc  =  (B„  X„  +  DTeuCd)  (1  -  Y2YxX„yl 

Kf  =  Y„cJ  +  BdDyd 

Kcl  =  -<7'2DydBd  Xtt  +  Cy)  (I  -  y2Y 

K„  =  Y'2Y„CI  Deu  +  Bu 

Xg,  and  Y^,  are  the  solutions  to  the  ARE’s 

(A  -  BuDduCe)TX(30  +  X„,(A  -  BuDduCe)  +  XM(7‘2BdBd  -BUB*  )XM 

+  Cd  (I  -  DeuDdu)T(l  -  DeuDTeu)Ce  =  0 

(A  -  BdDydCy) Y gg  +  Y„(A  -  BdDydCy)T  +  Yx(Y2c]  C.-C*  C,)Y„ 

+  Bd(I  -  DydDyd)(I  -  DTy„Dyd)TBj  =  0 

and  the  freedom  parameter  Q  is  given  by 

Q£*h„  llQlL  <  y 

There  are  three  conditions  that  must  be  met  in  order  for  this  parameterization 
to  be  valid: 

i)  Hx  6  dom(Ric)  with  X*  =  Ric(Hx)  >  0 

ii)  Hy  6  dom(Ric)  with  Y*  =  Ric(Hy)  >  0 

iii)  P(Y«XW)  <  72 
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where  Hx  and  Hy  are  the  associated  Hamiltonian  matrices  of  their  respective 
Riccati  equations,  and  p(  Y^XJ  is  the  spectral  radius  of  the  matrix  Y^X*.  If 
any  of  these  three  conditions  arc  not  met,  y  is  too  low  (that  is,  it  is  below  70). 
This  is  where  the  iterative  nature  of  HM  optimization  is  seen.  y0  must  be 
approached  iteratively,  and  7  can  be  made  to  be  arbitrarily  close  to  the  optimal 
value. 

Note  that  the  simplest  controller  in  this  parameterization  is  given  when  Q=0. 
This  is  known  as  the  central  compensator.  In  this  case,  the  compensator  is 
unique  and  its  order  is  equal  to  the  order  of  the  plant.  Other  higher  order 
compensators  are  in  the  set  of  solutions  (for  nonzero  values  of  Q),  but  the 
maximum  order  required  for  this  problem  is  the  order  of  the  plant. 

Without  rigorous  development,  consider  briefly  the  Ha  optimal  case.  When 
7  gets  close  to  yQ,  the  Ha  suboptimal  controller  (at  least  in  a  SISO  problem)  has 
a  singular  value  plot  that  is  completely  flat  from  low  frequency  out  to  a  high 
frequency  roll-off.  In  other  words,  it  looks  like  an  all  pass  filter  with  a  single 
high  frequency  pole.  As  7  approaches  70,  this  pole  moves  toward  infinity.  At 

optimal  (that  is,  7  =  70),  the  pole  is  at  infinity  and  the  controller  actually 
drops  rank  [Rid91b].  Therefore,  the  optimal  order  for  H*  optimal  compensators 
is  no  greater  than  full  order  minus  one. 
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2.4  Minimum  Entropy 

The  concept  of  the  entropy  of  a  control  system  is  very  non-intuitive  and  is 
not  directly  related  to  the  traditional  ideas  of  thermodynamic  randomness.  Let 
G  G  RLoo  and  jj  G  |j  <  7.  The  entropy  at  infinity  of  G(s)  is  defined  by 


s  lim 

' 

f  00  lnjdet[I  -  7  2  G(ju)  *  G(j«)]  | 

2 

So 

doj 

s0-°° 

2ir  J 

2  2 
so+w 

Note  from  examination  of  this  definition  that  the  entropy  has  the  following 
properties: 

i)  I[G(s),7]  *  0 

ii)  I[G(s),7]  =  0  iff  G(s)  =  0 

iii)  I[G(s),7]  <  00  iff  G(oo)  =  0  [MG90,8-1 1] 


While  this  definition  does  not  appear  to  have  much  practical  usefulness  in  this 
form,  entropy  is  a  useful  quantity  and  has  some  interesting  relationships  to  the 
2-norm  and  00 -norm.  In  particular,  entropy  is  an  upper  bound  to  the  2-norm  and 
is  equal  to  2-norm  as  7  gets  very  large.  That  is,  if  G  G  RH2  and  7  is  such  that 
II G  |  „  <  7,  then  I[G(s) ,7]  >  ||G(s)|||  and  I[G(s),«]  =  ||  G(s)  || |  <[MG90,12], 
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Theorem  2.4.4).  Also,  entropy  is  relatively  easy  to  compute  in  state  space 
despite  its  formidable  definition.  Let  G  G  RH2,  ||  G  ||  <»  <  7,  and 
G(s)  =  CCsI-A)'^.  Then 

I[G(s),7]  =  tr  [QCtC] 

where  Q  =  QT  >  0  is  the  stabilizing  solution  to  the  ARE 

0  =  AQ  +  QAt  +  7‘2QCtCQ  +  BBt  [MG90, 52-54]  Lemma  5.3.2 

Now,  consider  minimizing  the  closed-loop  entropy  of  a  system.  Make  all  the 
same  assumptions  that  were  made  for  the  H^,  suboptimal  parameterization  and 
let  7  >  70.  Then,  minimizing  the  closed-loop  entropy  I[Ted(s),7]  over  the  set 
of  all  admissible  compensators  K(s)  such  that  ||  Ted  ||  »  <  7  results  in  the  central 
H^,  controller  at  that  7  [MG90].  The  order  of  the  central  controller  has 
already  been  shown  to  be  no  greater  than  the  order  of  the  plant.  Therefore,  the 
minimum  entropy  problem  has  an  optimal  order  that  is  no  greater  than  the  order 
of  the  plant. 
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Table  2-1  gives  a  summary  of  the  optimal  orders  of  the  compensators  that  are 


solutions  to  the  optimization  problems  just  discussed. 


Table  2-1.  Optimal  Compensator  Orders 


Optimization  Optimal  Order  of 

Problem  Compensator 


H2  (opt) 

0 

(full-state  feedback) 

<  n 

(output  feedback) 

Ha  (sub-opt) 

<  n 

Ha,  (opt) 

<  n-; 

l 

I  [Ted,7l 

<  n 

note:  n  =  order  of  plant 
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III.  General  Development  of  the 
Mixed  FU/Hec  Problem 


3.1  Problem  Statement 

Consider  again  the  general  form  of  the  linear  time-invariant  feedback  control 
system  for  the  mixed  H2/Hoe  optimization  problem  as  introduced  in  Chapter  1, 
reproduced  in  Figure  3-1. 


Figure  3-1.  Mixed  System  Block  Diagram 


The  exogenous  input  d  is  assumed  to  be  a  deterministic  signal  of  unknown  but 
bounded  energy.  The  input  w  is  assumed  to  be  a  zero-mean  white  Gaussian  noise 
of  unit  intensity.  The  signals  e  and  z  are  controlled  outputs,  and  may  be  equal, 
dependent,  or  independent.  The  signals  u  and  y  are  the  control  inputs  and 
measured  outputs,  respectively.  All  weighting  functions  required  to  obtain  this 
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general  form  are  included  in  the  generalized  plant  transfer  function  P.  The  state 
space  equations  for  P  are  given  by: 


dr/dt  =  Ax  +  B  dd  +  Bwvv  +BU« 

(3.1a) 

e  =  +  Dedrf  +  Dewvv  4-  Deu« 

(3.1b) 

z  =  Cj  +  D  zdd  +  Dzwvv  +  Dzu« 

(3.1c) 

y  =  CyX  +  D ydd  -1-  Dywvv  +  Dyu« 

(3. Id) 

The  transfer  function  matrix  for  the  open-loop  plant  P  can  be  partitioned  as 


P  = 


* 

ed 

P 

*  ew 

Peu 

> 

zd 

P 

zw 

Pzu 

yd 

P 

yw 

Pyu 

Closed-loop  transfer  functions  from  w  to  z  and  from  d  to  e  are  given  by  the  lower 
linear  fractional  transformations 


T 


ZW 


PZW  +  Pzu  K  [I-Py.K]-1  Pyw 


Ted  =  Ped  +  Peu  K  [I-P^K]’1  Pyd 
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Now,  assume  the  plant  P  satisfies  the  following  conditions: 


(0 

Dzw  =  0 

(ii) 

E 

>.d  =  ° 

(iii) 

Dyu  =  0 

(iv) 

(A,BU)  stabilizable 

(v) 

DTz„Dz„, 

D  DT 

^yw^yw 

(vi) 

DeuDeu’ 

DydDyd 

(vii) 

A  —  jco  I 

Bu 

Cz 

Dzu  _ 

(viii) 

A-jodl 

K 

Dyw 

(ix) 

A-jodI 

Bu 

ce 

Deu 

(x) 

A-jodI 

Bd 

cy 

°yd 

(Cy,A)  detectable 
full  rank 
full  rank 

full  column  rank  for  all  u  €  R 

full  row  rank  for  all  w  €  R 

full  column  rank  for  all  a>  E  R 

full  row  rank  for  all  o  E  R 


The  rationale  for  these  assumptions  are: 

i)  Required  for  a  well-posed  H2  problem.  If  Dzw  ^  0,  the  2-norm  of 
Tzw  will  be  infinite,  regardless  of  K. 

ii)  Not  required  for  finite  oo-norm,  but  greatly  simplifies  the 
development. 


3-3 


iii)  Not  required,  but  greatly  simplifies  the  development  and  is 

representative  of  physically  realizable  plants. 

iv)  Necessary  for  the  existence  of  any  stabilizing  controller. 

v)  -  x)  Ensures  the  pure  H2  and  optimization  problems  each  have  an 

admissible  solution. 

The  following  definitions  will  be  used  throughout: 

To  s  inf  llTedlloo 

K  adm 

«o  s  inf  II  Tzw  II 2 

K  adm 

K2oPt  s  unique  K(s)  that  gives  ||TzvJ2  =  a0 
72  3  II  Ted  II  oo  when  K(s)  =  K2opt 

Kmix  s  a  K(s)  that  solves  the  mixed  H^H^,  problem  at  some  y 
y'  ■  |Ted|L  when  K(s)  -  Kmix 

«'  -  II TJW II 2  when  K(s)  =  Kmix 

n  as  order  of  the  plant 
nc  s  order  of  the  compensator 
n*  =  optimal  order  of  Kmix 
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Before  getting  too  deep  into  the  development  of  the  mixed  problem,  briefly 
survey  the  big  picture.  Consider  plotting  the  2-norm  of  Tzw  versus  the  oo-norm 
of  Ted.  Figure  3-2  shows  a  generic  plot  with  some  key  points  identified. 


Figure  3-2.  Generic  H2  versus  Plot 


This  plot  serves  to  illustrate  the  boundaries  of  the  mixed  problem.  Note  that  a0 
is  the  lowest  achievable  2-norm  by  any  compensator  since  it  comes  from  the 
unconstrained  H2  problem.  Likewise  y0  is  the  minimum  achievable  oo-norm. 
Therefore,  the  dashed  lines  represent  limits  of  achievable  H2  and  Hw 
performance  within  which  the  mixed  solutions  must  lie.  If  7  is  lowered  beyond 
70,  no  solution  will  exist.  Note  also  that  as  the  Ha  optimal  solution  is 
approached,  the  2-norm  becomes  infinite.  Actually,  there  are  possible  cases 
where  the  optimal  H,,,  controller  has  a  zero  state  space  D  matrix,  thus  allowing 
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a  finite  2-norm.  These  cases  are  rare  and  will  not  be  considered,  even  though 

the  development  to  follow  would  still  hold.  Finally,  as  7  is  increased  to  y2,  the 

solution  becomes  the  unique  H2  unconstrained  optimal  solution.  The  solution  will 

remain  the  same  unique  H2  optimal  solution  if  7  is  increased  further  because  the 

00 -norm  constraint  is  effectively  removed. 

Now,  the  general  mixed  H2/H00  optimization  synthesis  problem  is: 

Determine  the  set  of  admissible  compensators  K(s)  such  that 

inf  1  Tzw  1 2,  subject  to  the  constraint  ||  Ted  ||  <7 

K  adm 

is  achieved. 


First,  this  problem  statement  needs  to  be  turned  into  a  mathematical  statement 
that  can  be  manipulated.  Begin  with  a  definition  of  the  compensator.  A  state 
space  description  of  the  compensator  in  Figure  3-1  is 

dbcc/dt  =  +  B^y  (3.2a) 

u  =  Ccxc  +  Dcy  (3.2b) 


It  is  easily  shown  that  Dc  must  be  zero  in  order  to  have  a  finite  Tzw  [Rid91a,92], 
so  the  assumption  that  K(s)  is  strictly  proper  is  made  with  no  loss  of  generality. 


Using  the  control  law  u  =  K(s)y,  form  the  closed-loop  system  with  the 


augmented  state  vector 


by  combining  Equations  (3.1)  and  (3.2). 


3-6 


The  closed-loop  system  can  now  be  written  in  state  space  form  as 


—  =  AX  +  BAd  +  B 
dt  d 


W 


VV 


(3.3a) 


e  =  Cex  +  Dev,w 


(3.3b) 


Z  -  czl  -  Dzdd 


(3.3c) 


where  the  closed-loop  (tilde)  matrices  are  given  by 


A 

BuCc 

A  = 

BcCy 

Ac 

Bd 

Bw 

Bd  = 

BcDyd 

Bw  = 

BcDyw 

ce  =  [  c«  DeuCc  ]  cz  =  [  C,  DznCc  ] 


The  closed-loop  transfer  functions  from  w  to  z  and  from  d  to  e  can  now  be 
written  as 


Ted  =  Ce  (si  -  A)'1  Bd 
Tzw  =  Cz(sl  -  A)’1  Bw 
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Now  that  the  closed-loop  system  has  been  expressed  in  terms  of  the  unknown 
compensator  matrices  (Ac,  Bc,  Cc),  the  cost  objective,  stability  requirement,  and 
constraint  need  to  be  recast  as  mathematical  functions.  Begin  with  the  cost 
objective  and  stability  requirement.  It  is  a  standard  result  that  the  2-norm 
squared  of  a  transfer  function  can  be  found  as  a  function  of  the  solution  to  a 
Lyapunov  equation.  Specifically,  if 


!  A  !  B 


w 


Tzw(s)  - 


Cz  ,  0 


e  RH 


(which  is  true  with  the 
given  assumptions) 


then 

I TIW  ||  |  =  trlQ^CJ 

— 

where  Q2  =  Q2  ^  0  is  the  solution  to  the  Lyapunov  equation 

AQ2  +  Q2At  »  BWB«  =  0 

Note  from  Lyapunov  theory  that  this  symmetric,  positive  semidefinite  solution 
only  exists  when  A  is  stable  [Won85,283],  so  the  requirement  of  finding  a 
stabilizing  compensator  will  be  automatically  satisfied  if  the  solution  of  this 
equation  is  enforced. 


Now  consider  the  «-norm  constraint.  It  can  be  shown  that  an  <»-norm 


bound  on  a  transfer  function  can  be  guaranteed  by  requiring  the  solution  to  a 
Riccati  equation.  Specifically,  for 


and  assuming  A  is  stable  (which  is  enforced  by  the  Lyapunov  equation),  then  if 

T 

there  exists  a  Qw  =  Q*  >  0  that  satisfies  the  algebraic  Riccati  equation 


AQo.  *  Q„At  *  7~2Q»CeTCeQ.  .  gdBdT  =  0 


then 

IITJL  £  7 


[Rid91a, 55-57]  Thm  2.5.17 


Finally,  the  mixed  problem  can  now  be  restated  as: 

Determine  (Ac,  Bc,  Cc)  that  minimizes  the  cost  function 

J(AC,  Bc,  Cc)  =  tr  [Q?C^CZ]  (3.4) 

where  Q2  is  the  real,  symmetric,  positive  semidefinite  solution  to 

AQ2  +  Q2At  +  BW§J  =  0 
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and  such  that  (for  a  given  7) 


AQ»  *  Q„AT  ♦  7'2Q»CeTC,Q»  *  M7  -  0 

has  a  real,  symmetric,  positive  semidefinite  solution. 

Note  that  the  variables  in  the  problem  are  Ac,  Bc,  Cc,  Q2,  and  Q^.  All 
other  matrices  are  known  from  the  system  model,  and  7  is  given  in  the  problem 
statement. 

3.2  First-Order  Necessary  Conditions  for  the 
General  Mixed  Problem 

In  order  to  derive  the  necessary  conditions  for  the  mixed  I^/H,*,  optimization 
problem,  the  method  of  Lagrange  multipliers  is  used.  In  this  method,  a 
constrained  optimization  problem  is  transformed  into  an  unconstrained  problem 
by  adjoining  the  constraint  equations  to  the  original  cost  function  and  forming  a 
new  cost  function,  called  the  Lagrangian.  The  first  variations  of  this  Lagrangian 
are  then  set  equal  to  zero  to  find  the  necessary  conditions  for  a  minimum  (second 
variations  give  the  sufficient  conditions).  The  constraint  equations  (which  must 
be  of  the  form  f(x,y,...)  =  0)  are  adjoined  by  a  parameter  called  the  Lagrange 
multiplier.  This  Lagrange  multiplier  is  then  another  variable  which  must  be 
solved  for.  If,  during  the  solution  of  the  problem,  it  is  shown  that  the  Lagrange 
multiplier  must  be  zero,  it  physically  means  that  the  constraint  is  not  in  effect 
and  can  be  removed. 
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For  this  problem,  there  are  two  constraint  equations  to  be  adjoined.  Since 
they  are  matrix  equations,  the  Lagrange  multipliers  must  also  be  matrices. 
Define  these  two  matrices  in  partitioned  form  as 


X1 

*12 

Yi 

N> 

_ 1 

X  = 

1 

X 

Z  H 

X2 

Y  = 

yt 

1 12 

Y2_ 

These  matrices  do  not  have  to  be  defined  as  symmetric.  However,  they  will  be 
shown  to  be  real  symmetric  solutions  to  Lyapunov  equations,  so  for  simplicity 
they  will  be  defined  to  be  symmetric  here  (Xj,  X2,  Yt,  and  Y2  are  therefore 
symmetric). 

The  Lagrangian  can  now  be  formed  and  the  mixed  Hj/H^  Lagrange 

multiplier  problem  can  be  stated.  Assume  Q2  =  Q2  ^  0  and 
T 

Qoo  =  Qoo  -  0-  Minimize  the  Lagrangian 

2  =  tr[Q2c!cj  +  tr{[AQ2  +  Q2AT  +  BWB*]X} 

+  tr([AQ„  +  Q„At  +  r2Q„,cI  CeQ„  +  B„Bjm  (3.5) 

Note  that  if  Y  =  0,  the  oo-norm  bound  is  trivially  satisfied.  That  is,  the 
constraint  is  inactive.  If  the  solution  lies  on  the  boundary  of  this  constraint,  then 

Y  *  0.  For  the  special  case  7  =  y2,  the  solution  lies  on  the  boundary  and 

Y  ssO. 
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Now  set  the  first  variations  of  this  function  equal  to  zero.  That  is,  take  the 
partial  derivatives  of  %  with  respect  to  all  the  variables  (Ac,  Bc,  Cc,  Q2,  Qoo. 
X,  and  Y)  and  set  them  equal  to  zero.  These  are  derivatives  of  scalars  with 
respect  to  matrices  for  which  the  following  formulae  apply: 


a  tr(AXB)  _  tbt 

ax 

a  tr(AXTB)  _  BA 

ax 

9  tr(AaxBXC)  =  aTcTxT®t  +  btxtatct 

^ir..(AX.BxTc)  =  AtCtXBt  +  CAXB 

ax  [AS65] 


The  derivatives  with  respect  to  Q2,  Q,*,,  X,  and  Y  are  no  problem  because  they 
appear  explicitly  in  the  Lagrangian  in  its  present  form.  However,  the 
compensator  matrices  Ac,  Bc,  and  Cc  do  not  appear  explicitly,  so  the  Lagrangian 
must  be  multiplied  out  and  expanded  into  its  constituent  sub-blocks.  Define  the 
sub-blocks  of  Q2  and  QM  as  (note:  these  are  symmetric  matrices) 


q2 


Qi 

Ql2 

Qa 

Qab 

Q12 

Q2 

Qoo  - 

_QaTb 

Qb_ 
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The  off-diagonal  terms  in  these  partitions  are  not  square  unless  the  order  of  the 
compensator  is  equal  to  the  order  of  the  plant.  In  general,  the  sizes  of  the 
matrices  within  the  partitions  of  A,  Q2,  Qoo>  X,  and  Y  are 


(n  +  nc)  x  (n  +  nc) 


n  x  n  :  n  x  nc 

&c x  n  j  Be  1  nc 


In  addition,  to  reduce  excessive  notation,  make  the  following  simplifying 
definitions: 


Va  =  BdBd 

Vab  =  BdDTyd 

Vb  =  DydDTyd 

Ra  =clce 

R  =  C>eu 

Rb  ~  ^eu^eu 


Vj  =  bwb l 
V 12  =  BwDyw 
V2  =  DywDyw 

Ri  =  C*CZ 
r12  =  Cz  Dzu 
R2  "  DzuDzu 
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After  substituting  all  the  expressions  for  the  closed-loop  system,  Q2,  Q«>,  X,  Y, 
the  compensator,  and  the  simplifying  definitions,  the  expanded  Lagrangian  is 


X  =  tr{  Q,R,  +  Qi2cI  R^2  +  Q*2Ri2Cc  +  Q2C*  R2Cc  +  AQ,X, 

+  Q^X,  +  B.C.tfljX,  +  Q12Cc  bJ  X,  +  V,X,  +  AQ,2X^2  +  Q12Aa  X*2 
+  BuCcQ2X^2  +  Q,c,  Bb  x*2  +  V12B^xT2  +  BcC,Q,X12 
+  acqy2x12  +  q!2atx12  +  Q2cl  b*  xl2  +  bcvy2x12  +  acq2x2 

+  Q2AYX2  +  BcCyQ,2X2  +  X2  +  BcV2Bb  X2  +  AQ.Y, 

+  QaATY,  +  B„CcQtY,  +  Qatcl  Y,  +  V.Y,  +  AQabY*2  +  QabA*  Y*2 
+  BuCcQ„Y^2  +  QaCy  Bb  YT2  +  VabBa  Y\2  +  BcC,QaY12  +  AtQTabYl2 
+  QlbATY12  +  QbCa  Bj  Y12  +  BcVabY12  +  AcQbY2  +  QbAa  Y2 
+  BcCyQab  Y 2  +  Qa,Cy  BTc  Y2  +  Bc  V2Bb  Y2  +  y  \  QaRaQaY , 

+  QabCl  RT,bQ.v,  +  Q,RabccQab  y  i  +  Q«bc l  RbccQTabY,  +  QaRaQabY*2 
+  QabcI  RTabQabYa  +  Q.R,bCcQbYL  +  Qabcl  R„ccQbYY2 
+  QTabR»Q.Yi2  +  QbclRT,bQ,Y12  +  QTabRabccQTabY12 
+  Q„clRbCcQTabY12  +  QTabRaQabY2  +  Qbcl  R*„QabY2 
+  QTabRabCcQbY2  +  QbCb  RbCcQ„Y2)  }  (3.6) 
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Now,  the  partial  derivatives  of  the  Lagrangian  are  [Rid9 la, 97-98] 


=  X,T2Q12  *  X2Q2  ♦  Y,T2Qab  *  Y2Qb  =  0 


(3.7) 


-  X!T2Q,CyT  ♦  X2Q,T2CyT  +  X,T2V12  ♦  X2BcV2  *  Y,T2QaCyT 
*  Y2QaTbC,T  *  Y,T2Vab  *  Y2BcVb  =  0 


(3.8) 


*  B„TX,Q,2  *  B„TX12Q2  ♦  R,T2Qi2  *  R2C£Q2  *  B„TY,Qab 

*  BuTY12Qb  ♦  7‘2  [  RaTbQaY!Qab  ♦  RaTbQaY12Qb 

*  R,TbQ,bV,T2Qab  ♦  RaTbQ,bY2Qb  *  RbCcQtTbY,Qab 

+  RbCcQbYjT2Q>b  *  RbCcQaT„Yl2Qb  *  RbCcQ„Y2Qb  )  -  0 


d  % 
dx 

dj£ 

dQ2 

d  % 
dY 

d  % 

dQo. 


AQ2  +  Q2At  +  BWB^  =  0 
AtX  *  XA  +  CZTCZ  =  0 
AQ«  ♦  Q.A1  +  7'2QooCeTCeQoo  +  BdBdT  =  0 
[A  +  7_2QooCeTCe]T  Y  +  Y  [A  +  7'2QooCeTCe] 


(3.9) 

(3.10) 

(3.11) 

(3.12) 

0  (3.13) 
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These  are  a  set  of  seven  nonlinear,  coupled  matrix  equations  for  which  no 
analytical  solution  has  been  found.  Equations  (3.10)  and  (3.11)  are  Lyapunov 
equations  in  Q2  and  X,  respectively.  It  is  interesting  to  note  that  these  two 
equations  show  that  Q2  and  X  are  the  controllability  and  observability  gramians 
of  T2W.  Equation  (3.13)  is  a  Lyapunov  equation  in  Y,  but  it  has  no  constant 
term.  Therefore,  if  the  matrix  (A  4-  y  Q00CeCe)  is  stable,  the  only  solution 
to  (3.13)  is  Y  =  0  ([SZ70],  Thm  2.1).  Now,  Equation  (3.12)  is  a  Riccati 
equation  in  Q^,.  If  Q*  is  chosen  to  be  the  stabilizing  solution  to  this  Riccati 
equation  (which  is  typically  the  solution  that  is  wanted),  the  matrix  that  is 
stabilized  is  (A  +  7  QooCe  Ce).  This  means  that  the  Lagrange  multiplier  matrix 
Y  must  be  equal  to  zero,  which  means  that  the  00 -norm  constraint  is  not  active 
and  can  effectively  be  removed. 

This  apparent  contradiction  might  lead  one  to  question  the  initial  set-up  of 
the  problem.  However,  one  of  the  key  contributions  of  Ridgely’s  work  was  the 
recognition  that  the  neutrally  stabilizing  solution  (not  the  stabilizing  solution)  to 
(3.12)  is  the  solution  that  is  desired  in  order  to  obtain  an  on-boundary  mixed 
solution.  If  (A  +  7'2QoeCgCe)  is  neutrally  stable  (that  is,  it  has  at  least  one 
ju-axis  eigenvalue  with  the  remaining  eigenvalues  in  the  left  half  plane),  then 
there  are  an  infinite  number  of  non-zero  solutions  (Y)  of  varying  rank  ([SZ70], 
Lemmas  4. 1  and  4.2).  This  immediately  brings  up  the  question  of  uniqueness  of 
the  mixed  solution.  No  proofs  of  uniqueness  have  been  given,  and 

although  this  issue  was  touched  upon  in  this  research  while  trying  to  prove 
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optimal  order,  no  proofs  of  uniqueness  will  be  offered  here.  While  it  may  turn 
out  to  be  true  that  the  solution  is  indeed  unique,  it  would  not  be  surprising  to 
find  out  that  it  is  not,  since  both  the  H2-suboptimal  and  H^-suboptimal  (and 
optimal)  are  non-unique. 

Since  these  equations  are  intractable  analytically,  they  must  be  solved 
numerically  (Note:  even  though  closed-form  solutions  for  Ac,  Bc,  Cc  are  not 
available,  some  important  results  can  be  shown  from  these  equations.  Some  of 
the  key  theorems  that  describe  the  nature  of  the  mixed  problem  for  full  order 
controllers  are  given  later  in  Section  3.4).  There  are  software  programs  readily 
available  that  solve  Riccati  and  Lyapunov  equations.  However,  they  do  not 
return  neutrally  stabilizing  solutions  to  Riccati  equations  and  do  not  handle 
Lyapunov  equations  that  have  no  constant  terms.  The  equations  could  be 
numerically  handled  without  using  any  Riccati  or  Lyapunov  solvers,  but  the 
number  of  unknowns  becomes  unreasonable  very  quickly.  Also,  the  numerics 
become  erratic  near  the  solution  because  the  solution  lies  just  on  the  boundary 
of  non-existence  of  solutions.  Therefore,  even  to  solve  this  problem  numerically, 
modifications  are  needed.  This  is  the  motivation  for  setting  up  the  £u/>optimal 
mixed  ^/H^,  problem,  described  next. 
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3.3  First-Order  Necessary  Conditions  for  the 
Suboptimal  Mixed  Problem 

The  objective  now  is  to  develop  a  suboptimal  mixed  problem  that  has 

a  set  of  necessary  conditions  that  are  more  numerically  stable  and  that  approach 
the  optimal  solution.  Recall  the  statement  of  general  mixed  problem: 
Determine  (Ac,  Bc,  Cc)  that  minimizes  the  cost  function 
J(AC,  Bc,  Cc)  =  tr  [Q2C^Cz] 

where  Q2  is  the  real,  symmetric,  positive  semidefinite  solution  to 
AQ2  ♦  Q2At  +  BWB l  =  0 
and  such  that 

AQo.  +  Q„AT  *  7‘2Q»CeTC,Q„  *  B„BdT  =  0 
has  a  real,  symmetric,  positive  semidefinite  solution. 

Now,  for  the  suboptimal  mixed  problem,  define  a  new  cost  objective  (the 
constraint  equations  remain  unchanged): 

Jm(Ac,  Bc,  Cc)  =  (1-m)  tr[Q2Cz  Cz]  +  M  trlQ^C*  Ce]  (3.14) 

where 

n  E  R,  n  G  [0,1] 
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Notice  that  when  n  =  0,  the  problem  reduces  to  the  original  optimal  mixed 

— T  — 

problem.  The  additional  term,  Ce)  in  the  cost  function  is  partially 

variable  as  (CeCe)  could  actually  be  any  real,  symmetric,  positive  semidefinite 

~T  — 

matrix.  The  reason  for  choosing  (CeCe)  is  that,  for  n  ^  0,  the  term 
tr[QooCgCe]  is  the  entropy  of  Ted  ([Rid9 1  a,  1 1 8] ,  Thm  5.1.1).  This  is  very 
convenient  because  it  provides  an  excellent  starting  point  for  the  numerical 
algorithm.  When  fx  =  1,  the  problem  reduces  to  pure  entropy  minimization  for 
which  a  closed-form  solution  is  available  (it  is  the  central  HM  controller  at  the 
given  level  of  7). 

The  reason  for  adding  the  extra  term  into  the  cost  function  was  to  try  to 
modify  the  necessary  conditions  so  that  the  Lyapunov  equation  in  the  necessary 
conditions,  Equation  (3.13),  would  have  a  small  symmetric  positive  semidefinite 
constant  term.  This  would  enable  the  use  of  a  Lyapunov  solver  and  would 
require  the  term  (A  +  7  Q<»Ce  Ce)  to  be  stable  without  Y  =  0.  Hence,  the 
stabilizing  solution  to  Equation  (3.12)  would  be  desired  and  a  Riccati  solver 
could  also  be  used.  The  additional  term  does  indeed  accomplish  this. 

With  the  new  cost  function  defined,  the  suboptimal  mixed  Lagrangian 

becomes 

=  d-M)  tr[Q2C^CJ  +  M  trtQ^CJ 
+  tr{[AQ2  +  Q2At  +  B  wBj]X} 

+  trUAQoe  +  QoeAT  +  7'2Q»cICeQ00  +  BdBd  ]Y}  (3.15) 
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Now,  the  partial  derivatives  of  this  new  Lagrangian  are  [Rid91a,  116-117] 


=  X,T2Ql2  *  x2q2  +  Y,T2QaB  *  Y2Qb  =  0 


(3.16) 


a  % 


n  _ 


,T 


aB„ 


=  x^QjC;  ♦  x2q,‘2c;  +  x12v12  +  x2bcv2  +  Y12qac. 


♦  Y2QaTbCyT  *  Y,^Vab  *  Y2BcVb  =  0 


(3.17) 


4^  =  B„TX,Q12  *  BuTX12Q2  *  (1-^)R,T2Q12  »  (1-^)R2CcQ2 
dCc 

*  B„TY,Qab  <  BuTY12Qb  *  MRaTbQab  *  MRbCcQb 

*  r1  l  RaTbQaY,Qab  .  R.T„QaY12Qb  *  RaTbQabY,T2Qab 
+  R£OiQabY2Qb  +  Rb^cQabYlQab  +  RbCcQbYi2Qab 

*  RbCcQaTbY|2Qb  *  RbCcQbY2Q„  ]  =  0 


_ _ M 

ax 
a  % 

_ m 

aQ2 

a^ 

ay 

^Qoo 


aq2  +  q2at  +  b  wbI  =  o 

AtX  +  XA  +  (1  -t'.)CzTCz  =  0 

AQoe  -  Q»  AT  -  7*2QccCeTCeQoo  +  BdBdT  =  0 
[A+7-2Q06CeTCe]TY  +  Y  [A  +7_2QoecJ Ce]  +  MCeTCe 


(3.19) 


(3.20) 


(3.21) 


0 

(3.22) 
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It  is  immediately  obvious  that  the  clever  modification  of  the  cost  function  is 
very  helpful  for  setting  up  a  more  well-behaved  problem  for  numerical  solution. 
Equation  (3.22)  is  a  typical  Lyapunov  equation,  and  Equation  (3.21)  now  requires 
a  stabilizing  solution.  At  /x  =  0,  the  necessary  conditions  are  exactly  the  same 
as  the  ones  derived  for  the  optimal  mixed  problem.  At  ;x  =  1,  the  problem  is  the 
minimum  entropy  problem  for  which  a  closed-form  solution  can  be  found.  Now, 
all  that  is  required  in  the  numerical  solution  is  to  start  with  the  central 
controller  (with  /x  =  1)  and  successively  reduce  n  until  the  solution  converges  to 
the  optimal  mixed  solution  (which  it  does,  as  will  be  shown  in  the  next  section). 

3.4  Summary  of  Full  Order  Results 

The  suboptimal  mixed  problem  is  set  up  purely  for  the  sake  of  performing 
the  numerical  solution.  Even  though  the  optimal  mixed  problem  could  not  be 
solved  in  closed-form,  some  important  characteristics  of  the  solutions  have  been 
proven  for  compensators  of  equal  order  to  the  plant.  The  key  proofs  that  relate 
directly  to  the  increased  order  problem  are  summarized  here. 

Notice  that  until  now,  the  order  of  the  compensator  K(s)  has  not  been 
specified,  so  all  the  necessary  conditions  are  valid  for  any  order.  However,  as 
compensator  order  changes,  the  sizes  of  the  equations  change.  This  means  that 
in  order  to  continue  with  the  analysis  (either  numerically  or  analytically)  the 
order  of  the  compensator  must  be  chosen.  In  Ridgely’s  work,  the  order  was 
selected  to  be  full  order  (that  is,  nc  =  n).  The  majority  of  the  proofs  and 
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examples  given  in  this  section,  therefore,  are  given  under  the  assumption  of  a  full 
order  compensator. 

Theorem  3.4,1:  For  nc  of  any  order  and  7  <  70,  there  is  no  solution  to  the 
mixed  problem. 

Proof:  See  [Rid91a,  102],  Thm  4.2.2. 

Theorem  3.4.2:  Assume  nc  >  n,  and  7  is  selected  such  that  7  >  y2.  Then 

^mix  —  ^2opt 

ii)  a*  =  a0 

iii)  7*  =  72 

Proof:  See  [Rid91a,101],  Thm  4.2.1. 

Theorem  3.4.3:  Assume  nc  =  n,  and  7  is  selected  such  that  70  <  7  <  y2. 
Then  the  solution  to  the  mixed  Hj/H^  problem  lies  on  the  boundary  of  the 
00 -norm  constraint.  That  is,  Kraix  is  such  that  7*  =  7. 

Proof:  See  [Rid91a,104],  Thm  4.2.3. 

Theorem  3.4.4:  Assume  nc  =  n.  A  plot  of  a*  versus  7*  for  7  >  yQ  is 
monotonically  decreasing  with  7. 

Proof:  See  [Rid9 la,  123],  Thm  5.1.5. 
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There  are  also  some  important  theorems  for  the  suboptimal  mixed  problem. 
One  of  the  most  important  insures  that  the  suboptimal  problem  converges  to  the 
optimal  solution. 

Theorem  3.4.5:  Assume  nc  =  n,  and  7  >  70.  Then  given  by  (3.14) 
converges  to  the  optimal  mixed  Hj/H^  problem  as  ^  -*  0. 

Proof:  See  [Rid91a,  120],  Thm  5.1.2. 

One  the  best  ways  to  visualize  all  of  the  above  theorems  is  by  examining  a 
typical  plot  of  ||  Ted  ||  «  versus  ||Tzw||2. 


Figure  3-3.  Example  Full  Order  Mixed  Plot 


This  plot  is  taken  from  actual  data  from  Ridgely’s  full  order  SISO  example 
[Rid91a,  167].  Several  characteristics  are  immediately  evident.  All  solutions  lie 
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inside  of  (or  on)  the  a0  and  y0  boundaries.  The  2-norm  gets  very  large  near  yQ. 
The  mixed  curve  is  monotonically  decreasing  and  does  not  reach  aQ  until  y2. 
The  mixed  curve  literally  stops  at  the  (a0,y2)  point.  Also  shown  on  the  plot  are 
the  central  solutions  (Q  =  0).  Notice  the  large  improvement  of  the  mixed 
solutions  over  the  Q  =  0  solutions.  This  plot  is  very  encouraging  from  a 
designer’s  standpoint.  It  shows  that  a  significant  reduction  in  the  00-norm  can 
be  achieved  by  sacrificing  only  a  small  amount  of  H2  optimality.  Physically,  this 
means  that  there  are  regions  where  the  system  can  be  made  much  more  robust 
without  giving  up  very  much  H2  performance. 

Recall  that  all  the  above  results  assume  a  full  order  compensator.  What 
would  the  mixed  plot  look  like  for  higher  order  controllers?  Do  the  higher  order 
controllers  exhibit  the  same  characteristics?  Is  it  possible  to  increase  the 
compensator  order  enough  to  achieve  aQ  at  level  of  y  below  y2l  These  are  issues 
that  are  addressed  in  the  next  two  chapters. 
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IV.  Increased  Order  Compensators  in 
Mixed  Optimization 


4.1  Optimal  Mixed  Problem 

In  Chapter  III  the  general  mixed  I^/H^  problem  was  developed.  During  the 
development,  the  form  of  the  compensator  was  assumed,  but  its  order  was  not 
specified.  Then,  after  finding  the  necessary  conditions  for  a  minimum,  there 
came  a  point  where  the  order  of  the  compensator  had  to  be  fixed  before  the 
analysis  could  continue.  Results  for  full  order  compensators  have  been 
demonstrated;  however,  increased  order  controllers  have  not  yet  been  thoroughly 
addressed.  Due  to  the  dynamic  nature  of  solutions  under  the  assumption  of  state 
availability,  it  was  informally  conjectured  by  Rotea  &  Khargonekar  that 
compensators  of  order  greater  than  the  plant  would  be  required  in  the  output 
feedback  case  [RK91,310]  (this  conjecture  will  be  shown  to  be  true  for  certain 
choices  of  7).  The  issue  of  higher  order  compensators  is,  therefore,  well 
motivated.  This  chapter  extends  some  of  the  full  order  results  to  compensators 
of  order  greater  than  the  plant  for  the  output  feedback  case,  and  takes  some 
preliminary  steps  toward  answering  the  question  of  optimal  compensator  order 
for  the  mixed  problem. 
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First,  when  examining  the  mixed  problem,  it  will  be  helpful  to  define  three 
regions  based  on  the  chosen  level  of  7: 

Region  1:  0  <  7  <  70 

Region  2:  70  <  7  <  12 
Region  3:  y2  -  7  —  00 

As  seen  already  in  the  full  order  case,  the  nature  of  the  solution  is  highly 
dependent  on  7.  Therefore,  each  region  will  be  addressed  separately. 

4.1.1  Region  1.  Region  1  can  be  dismissed  immediately,  as  there  is  no 
controller  of  any  order  that  will  meet  the  00 -norm  constraint  since  the  set  of  all 
H^,  suboptimal  controllers  is  then  empty  (see  also  Theorem  3.4.1). 

4.1.2  Region  2.  Region  2  (70  <  7  <  72)  is  the  region  of  greatest  interest 
since  it  is  here  that  there  may  be  competing  objectives.  Recall  that  in  the  state 
feedback  case,  if  7  is  chosen  such  that  a  solution  exists,  the  absolute  minimum 
of  II  Tzw  II 2  could  be  achieved  while  meeting  the  |j  Ted  )j  ^  constraint.  Under 
output  feedback,  this  is  not  true  for  7  levels  lower  than  y2- 

Lemma  4,1.1:  Assume  nc  >  n,  and  7  is  selected  such  that  yQ  <  7  <  y2. 
Under  output  feedback,  it  is  not  possible  for  the  mixed  l^/H*  controller  to 
achieve  the  absolute  minimum  2-norm  of  Tzw.  That  is,  a*  ^  a0. 
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Proof.  Assume  a*  =  a0.  Then,  by  the  parameterization  of  all  H2-(sub)optimal 
compensators,  Q=0,  and  the  compensator  is  unique  (K2opt).  7*  is  identically 
equal  to  y2-  This  is  a  contradiction  for  7  <  y2.  Therefore,  if  7  <  y2,  then 
a*Va0  regardless  of  compensator  order.  ■ 

It  has  been  shown  that  for  full  order  compensators  (that  is,  nc  =  n),  the 
solution  to  the  mixed  problem  lies  on  the  boundary  of  the  00 -norm 

constraint.  In  other  words,  in  Region  2,  7*  =  7.  This  result  is  also  true  for 
higher  order  compensators. 

Lemma  4.1.2:  Assume  nc  >  n,  and  7  is  selected  such  that  yQ  <  y  <  y2.  The 
solution  to  the  mixed  Hj/H^  problem  lies  on  the  boundary  of  the  00 -norm 
constraint  (7*  =  7). 

Proof:  For  nc  =  n,  see  Theorem  3.4.3. 

For  nc  >  n,  begin  by  assuming  that:  a)  the  solution  is  off  the  boundary 
(Y=0),  and  that  b)  K(s)  is  a  minimal  realization.  With  the  assumption  that  the 
solution  is  off  the  boundary,  the  first  order  necessary  conditions  (with  the  sub¬ 
blocks  expanded)  are  given  by: 
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^12  Q 12  "*■  X2Q2  —  0  (4.1) 

xLQjCJ  +  X2q[2cJ  +  xT2V12  +  X2BcV2  =  0  (4.2) 

X{Ql2  +  X12Q2  +  r[2Q12  +  R2CcQ2  =  0  (4.3) 

AQ,  +  Q^7  +  BuCcqT2  +  Q12cI  +  Vt  =  0  (4.4) 

AQi2  +  Q12aT  +  BuCcQ2  +  Q.C]  BTC  +  V12bJ  =  0  (4.5) 

AcQ2  +  Q2Al  +  BcCyQl2  +  o[2Cy  BTc  +  BcV2bI  =  0  (4.6) 

AtX!  +  XjA  +  Cy  bJ  XTn  +  X12BcCy  +  Rj  =  0  (4.7) 

ATX12  +  X12Ac  +  Cy  BTc  X2  +  X!BuCc  +  R12Cc  =  0  (4.8) 

a!  x2  +  X2Ac  +  cj  B^  X12  +  x[2BuCc  +  CTc  R2Cc  =  0  (4.9) 


Note  that  these  equations  are  valid  and  must  be  satisfied  regardless  of  the 
compensator  order  that  is  chosen. 

Consider  equation  (4.6).  Since  Q2  >  0,  from  [KJ72, 147-148]  it  follows  that 

Q2  >  0  (4.10) 

Ql2  =  Ql2Q2  +  Q2 
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where  Q2+  is  the  Moore-Penrose  pseudoinverse  of  Q2.  Therefore,  (4.6)  can  be 
rewritten  as 

(Ac  +  BcCyQ12Q2+)Q2  +  Q2(Ac  +  BcCyQ12Q2+)T  +  BcDywDjwBj  =  0  (4.11) 

Now,  since  (AC,BC)  was  assumed  controllable  (minimal  realization),  by  ([Won85], 
Lemma  2.1)  it  follows  that  (Ac4-BcCyQi2Q2+,Bc)  is  also  controllable.  Further, 
since  Dyw  has  full  row  rank,  (Ac  +  BcCyQ12Q2+,BcDyw)  is  controllable.  Now, 
using  the  dual  of  ([Won85],  Lemma  12.2),  equations  (4.10)  and  (4.11)  imply 
Q2>0.  Therefore,  Q2_1  exists.  A  similar  argument  applied  to  (4.9)  implies  that 
X2'!  exists. 

Now,  examine  equation  (4.1),  rewritten  as 

X2Q2  =  'X12Q12  (4.12) 

For  the  case  of  nc  >  n,  X12  and  Ql2  are  non-square,  X12  E  Rnxnc  and 
Q12ERnxnc.  Therefore,  the  highest  rank  that  the  product  X12Q12  can  have  is  n, 
and  by  (4.12)  this  implies  rank(X2Q2)  <  n.  Therefore,  X2'1  (X2ERncxnc)  and/or 
Q2->  (Q2€R"c“  c)  do  not  exist.  However,  it  was  already  shown  that  both  X2'1 
and  Q2'1  exist;  thus,  a  contradiction.  At  least  one  of  the  assumptions  was 
violated.  There  are  three  possibilities,  all  of  which  yield  the  same  conclusion: 

1)  If  K(s)  is  a  minimal  realization,  then  Y  *  0.  Therefore,  the  solution  lies 
on  the  boundary. 
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2)  If  Y  *  0,  then  K(s)  is  not  a  minimal  realization.  It  immediately  follows 
(see  [Rid91a,106])  that  K(s)  is  just  a  state  space  transformation  of  the  unique  nth 
order  compensator  K2opt,  plus  arbitrary  pole-zero  cancellations.  Again,  the 
solution  lies  on  the  boundary. 

3)  If  both  assumptions  were  incorrect,  then  K(s)  is  not  a  minimal  realization 

and  Y  ^  0.  The  solution  still  lies  on  the  boundary.  ■ 

With  these  results  in  hand,  the  characteristics  of  Region  2  can  be  summarized 
as  follows: 

Theorem  4.1.1:  Assume  nc  >  n,  and  y  is  selected  such  that  yQ  <  y  <  y2- 
Then, 

0  “o  <  ^ 

(where  a*m  is  a*  for  an  mlh  order  controller) 
ii)  7*  =  7 
iii)  n*  >  n 

Proof:  i)  aQ  <  a*n  follows  immediately  from  Lemma  4.1.1.  If  the  higher 
order  Kmix  is  simply  a  non-minimal  realization  of  the  full  order  K^,  a*n  =a*n. 
However,  there  is  an  infinite  set  of  H2-suboptimal  controllers  for  a  >  aQ,  so 
further  reduction  of  a*  may  be  possible.  This  reduction  is  indeed  possible,  as 
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will  be  shown  by  example.  It  is  conjectured  that  condition  i)  can  be  strengthened 
to  a*n  <  a*n  in  general,  but  this  still  requires  formal  proof. 

“c 

Condition  ii)  was  shown  in  Lemma  4.1.2.  The  solution  lies  on  the  boundary 
of  the  oo-norm  constraint  for  increased  order  controllers. 

Condition  iii)  follows  from  i).  However,  to  prove  it  more  formally,  assume 
the  optimal  order  of  Kmix  is  n.  Then,  there  can  be  no  other  compensator  of  any 
order  that  can  achieve  a  value  of  ]j  Tzw  ||  2  lower  than  a*n.  However,  as  will  be 
shown  in  the  examples,  higher  order  compensators  can  be  used  to  achieve  a  lower 
a*  than  that  achieved  with  an  nth  order  controller.  This  contradicts  the 
assumption  that  the  optimal  order  must  be  n.  Therefore,  in  general,  the  optimal 
order  may  be  greater  than  the  order  of  the  plant.  In  fact,  if  the  conjecture  that 

a*n  <  a*n  is  shown  to  be  true,  iii)  can  immediately  be  strengthened  to  n*>n. 

■ 

4.1.3  Region  3.  For  Region  3  (y2  ^  y  —  the  oo-norm  constraint  is 
inactive  and  the  problem  reduces  to  an  uncons  J  H2  optimization  problem 
for  which  optimal  order  is  known. 

Theorem  4.1.2:  Assume  nc  >  n  and  y  is  selected  such  that  y  >  y2-  Then 

^mix  ~  ^2opt 

ii)  a*  =  0to 

iii)  y*  =  y2 

iv)  n*  =  n 
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Proof:  For  nc  =  n,  the  H 2  optimal  controller  is  unique  and  is  given  by  the 
results  in  Section  2.1  with  Q  =  0.  Using  this  compensator,  by  definition 
||  Tzw  ||  2  =  cx0  and  ||  Ted  ||  ^  =  y2.  Since  a0  is  the  absolute  minimum  achievable 
2-norm  of  Tzw,  and  since  the  constraint  ||  Ted  ||  ^  <  y  is  trivially  satisfied  when 
y  >  y2,  condition  i)  is  true  for  nc  =  n. 

For  nc  >  n,  K2opt  makes  ||TZW||2  =  aQ,  which  is  the  minimum  achievable 
value  of  ||  Tzw  ||  2  using  any  compensator  of  any  order.  Thus  it  must  be  the 
solution  and  condition  i)  is  true  for  nc  >  n. 

Conditions  ii)  and  iii)  follow  immediately  from  the  definitions  given  in 
Section  3.1. 

Condition  iv)  holds  because  K2opt  is  a  unique  compensator  of  order  n.  If  nc 
is  chosen  such  that  nc  >  n,  Kmix  must  be  nothing  more  than  a  state  space 
transformation  of  the  unique  compensator  plus  arbitrary  pole-zero 
cancellations.  ■ 

Some  extensions  from  the  full  order  case  to  the  increased  order  case  have 
now  been  made.  However,  while  the  optimal  order  of  the  mixed  solution  is  not, 
in  general,  the  order  of  the  plant,  the  ultimate  question  of  what  the  optimal  order 
is  remains  to  be  proven. 
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4.2  Suboptimal  Mixed  Problem 


As  described  in  Section  3.3,  due  to  difficulties  with  the  numerics  of  the 
necessary  conditions  for  the  optimal  mixed  problem,  a  suboptimal  approach  is 
used  when  solving  the  problem  numerically.  However,  since  this  is  different 
problem,  the  nature  of  the  suboptimal  problem  needs  to  be  addressed.  A  main 
concern  is:  as  ^  -*  0,  does  the  suboptimal  problem  converge  to  the  optimal? 
Simply  because  the  n  =0  case  reduces  to  the  optimal  mixed  problem  does  not 
imply  that  the  function  (Equation  (3. 14))  converges  to  the  function  J  (Equation 
(3.4))  as  n  approaches  zero.  It  is  possible  for  a  function  to  approach  a  certain 
value  in  the  limit  and  have  a  discontinuity  at  that  point  (that  is,  the  value  of  the 
function  at  the  point  is  completely  different  than  the  value  of  the  limit).  For  full 
order  compensators,  it  has  been  shown  that  the  suboptimal  problem  does  approach 
the  optimal  as  (see  Theorem  3.4.5).  However,  before  numerical  solutions 
are  obtained  for  higher  order  compensators,  this  same  convergence  needs  to  be 
demonstrated  for  nc  >  n. 

Lemma  4.2,1:  Let  G(s)  =  C(sI-A)'1B  with  A  stable.  Define  ||  G(s)  [|  ^  s  7n  and 
let  7  >  7n.  Define  s  =  y  -  yn  >  0.  Then  as  e  -*  0,  the  value  of  the  entropy, 
I[G(s),7]  (which  is  singular  at  e  =  0),  converges  to 
tr[XnCTC] 

where  Xn  is  the  neutrally  stabilizing  solution  to  the  ARE 
AXn  +  XnAT  +  7n-2XnCTCXn  +  BBt  =  0 
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Proof.  Immediate,  frc.n  the  fact  that  Theorem  5.1.1  and  its  proof,  given  in 
[Rid9 la,  1 18-1 19],  are  completely  general.  Even  though  the  assumption  that 
nc  =  n  was  made  prior  to  this  theorem,  the  theorem  as  stated  is  not  a  function  of 
compensator  order.  ■ 

This  result  is  very  important  for  the  suboptimal  mixed  problem.  It  is  well  known 
that  the  entropy  at  £  =  0  is  infinite.  However,  until  this  theorem  was  given,  it 
was  generally  accepted  in  the  literature  that  the  entropy  continuously  approaches 
infinity  in  the  limit  as  e  -*  0.  If  this  was  true,  then  the  added  term  trfQ^Cg  Ce] 
in  the  suboptimal  mixed  cost  function  would  become  infinite  as  y. 

approaches  zero,  because  the  neutrally  stabilizing  solution  to  the  related  Riccati 
equation  (3.12)  is  required  in  the  optimal  mixed  solution.  However,  since  the 
entropy  approaches  a  finite  value  in  the  limit,  it  is  possible  for  the  suboptimal 
mixed  problem  to  converge. 

Theorem  4.2.1:  Assume  nc  >  n,  and  7  >  70.  Then  given  by  (3.14) 
converges  to  the  optimal  mixed  H2/H00  problem  as  n  -*  0. 

Proof  For  nc  =  n,  see  Theorem  3.4.5. 

For  nc  >  n,  proof  is  immediate  upon  recognition  that  Theorem  5.1.2  and  its 
proof  given  in  [x< id9 1  a,  1 20]  are  completely  genera!  and  do  not  depend  on 
compensator  order.  The  transfer  function  discussed  in  the  proof  is  the  closed- 
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loop  transfer  function  Ted,  and  no  restrictions  are  placed  on  its  order.  For 
convenience,  a  rough  outline  of  the  proof  is  as  follows: 

-  Ted  is  stable  and  strictly  proper  at  all  values  of  n 

-  it  follows  that  the  second  term  in  (3.14)  is  I[Ted,7]  for  ^  ^  0 

-  from  Lemma  4.2.1  the  entropy  converges  to  a  finite  bound  as  /x  -*  0 

-  the  term  n  MQ^CgCJ  therefore  approaches  zero  as  n  -*  0 

-  thus,  JM  -*  J  as  n  -*  0  ■ 

Even  though  no  rigorous  proofs  were  required  to  extend  the  full  order  case 
to  the  higher  order  case  in  proving  the  convergence  of  the  suboptimal  mixed 
problem,  this  was  an  important  issue  that  needed  to  be  shown.  The  suboptimal 
mixed  problem  does  indeed  converge  to  the  optimal  for  higher  order 
compensators.  Having  shown  some  key  theoretical  results,  now  consider  the 
numerical  solution  of  the  mixed  H2/Hcx)  problem. 
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V.  Numerical  Solution 


5.1  Davidon-Fletcher-Powell  Algorithm 

Since  a  closed-form  solution  of  the  mixed  problem  is  not  available, 

the  problem  must  be  solved  numerically.  There  are  many  numerical  optimization 
algorithms  available.  Ideally,  a  second  order  gradient  method  (Newton’s 
procedure)  might  be  desired,  but  by  definition,  this  type  of  algorithm  requires  the 
second  partials  of  the  function  being  minimized.  For  the  mixed  H2/H(X  problem, 
these  second  derivatives  are  fourth-order  tensors.  Therefore,  due  to  memory 
limitations  and  excessive  execution  times  (as  discovered  by  [Rid91a]),  this 
approach  was  avoided.  Instead,  an  algorithm  developed  by  Davidon  and  further 
described  and  refined  by  Fletcher  and  Powell  was  used  (see  [Fox71 , 104-109]). 
The  Davidon-Fletcher-Powell  (DFP)  algorithm  is  a  very  powerful  quadratically 
convergent  first-order  method.  It  does  not  require  the  second  derivatives,  but  it 
does  calculate  estimates  of  these  derivatives.  These  estimates  are  used  to  form 
a  variable  metric  which  is  improved  with  each  iteration.  The  flow  diagram 
shown  in  Figure  5-1  outlines  the  basic  progression  of  the  algorithm.  X  is  a 
column  vector  containing  the  unknown  variables.  For  the  mixed  problem, 

these  are  the  Ac,  Bc,  and  Cc  matrices  stretched  into  a  vector.  F(X)  is  the 
function  being  minimized.  For  the  mixed  problem,  F  s  JM  .  The  gradient, 
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Figure  5-1.  DFP  Flow  Diagram 

7F(X),  for  the  mixed  problem  are  the  partial  derivatives  of  the  Lagrangian  3?^ 
with  respect  to  Ac,  Bc,  and  Cc  stretched  into  a  column  vector.  S  is  the  matrix 
that  specifies  the  direction  that  the  current  guess  of  X  needs  to  move  in  order  to 
reduce  F(X).  Typically,  in  a  true  second-order  method,  S  =  -  J'1  VF(X),  where 
J  is  the  Hessian  (second  derivative  matrix).  As  already  mentioned,  DFP 
estimates  J'1.  This  estimate  is  a  symmetric,  positive  definite  matrix  defined  as 
H.  Thus,  in  DFP,  S  =  -  H  VF(X).  k  is  the  metric  that  controls  the  size  of  the 
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step  tc  be  taken  in  the  S  direction.  M  and  N  are  matrices  that  are  used  to 
calculate  H,  and  will  be  defined  shortly. 

The  iteration  proceeds  as  follows: 

1)  Begin  with  an  initial  X  vector  and  H.  The  initial  X  that  was  used  was 
the  initial  guess  of  the  compensator.  Section  5.2  discusses  how  this  initial  guess 
was  determined.  The  initial  H  was  chosen  to  be  the  identity  matrix;  thus  the 
initial  direction  of  movement  was  simply  the  negative  of  the  gradient  at  that 
point. 

2)  Compute  Xq+1  =  Xq  +  x*qSq,  where  /c*q  minimizes  F(Xq  +  *qSq).  The 
method  used  for  finding  k*  was  a  simple  one-dimensional  search  involving  a 
doubling/bisection  technique.  For  an  outline  of  this  technique  see  [Rid9 1,1 30] . 
More  sophisticated  methods  were  attempted  in  order  to  speed  up  execution  times 
(this  is  where  most  of  the  computations  are  required);  however,  this  method 
proved  to  be  the  most  reliable  in  the  face  of  such  a  difficult  function. 

3)  If  the  algorithm  has  not  converged,  compute  Hq+1  =  Hq  +  Mq  +  Nq, 
where 

Y q  =  VF(Xq+1)  -  VF(Xq) 


M, 


S  ST 

q  t 
S  1  Y 

q 


N, 


(H,Y„)(HqY,)T 

YJH„Y, 
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4)  Finally,  calculate  the  new  S  by  Sq+1  =  -  Hq  +  1  VF(Xq+1)  and  return  to 
the  k*  calculation  step  with  the  new  X  and  S. 

Convergence  of  the  solution  is  declared  when 


VF,THqVFq 
I  F(X_)  | 


<  e 


[Fox71, 104-109] 


where  e  is  a  small  positive  number  that  is  defined  by  the  user  prior  to  execution 
(values  of  10'6  were  typical). 

This  algorithm  was  coded  into  FORTRAN  to  run  on  a  VAX  mainframe.  The 
FORTRAN  program  is  given  in  the  appendix.  Note  that  in  the  code,  the  variable 
k  is  called  a.  This  is  because  the  program  was  written  to  be  consistent  with  the 
nomenclature  given  in  [Fox71].  However,  in  order  to  avoid  confusion  with  the 
definitions  of  a  given  throughout  this  work,  the  discussion  in  this  section  uses 
the  nomenclature  k.  A  PRO-MATLAB™  version  of  DFP  was  used  by  Ridgely  in 
his  work.  The  FORTRAN  version  was  found  to  be  about  7-10  times  faster. 


An  outline  of  the  process  for  using  the  DFP  program  to  solve  the  mixed 
^/H,*  problem  is  as  follows: 

1)  Select  the  desired  compensator  order  and  level  of  7  and  determine  an 
initial  guess  for  the  mixed  compensator.  As  will  be  shown  in  the  following 
section,  tnis  was  not  a  trivial  task.  The  program  is  extremely  sensitive  to  the 
starting  point. 
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2)  Start  with  a  relatively  high  value  for  n  (usually  between  0.1  and  0.5)  and 
run  the  program  until  it  converges. 

3)  Use  the  compensator  found  in  step  2)  as  a  new  initial  guess  and  re-start 
the  program  at  a  lower  value  of  /x. 

4)  Continue  the  process  of  reducing  n  until  the  solution  converges 
acceptably  close  to  the  optimal  mixed  solution  (values  of  around  10~4  were 
typically  attained  for  p.).  Three  items  were  checked  when  determining  if  the 
solution  was  "close  enough."  First,  the  2-norm  of  Tzw  would  be  checked  to  see 
if  it  was  still  becoming  smaller.  Generally,  when  convergence  was  declared, 
reductions  could  only  be  seen  in  the  fifth  or  six  decimal  place.  Second,  since  it 
was  known  that  the  solution  has  to  lie  on  the  boundary  of  the  c»-norm  constraint, 
the  c» -norm  of  Ted  was  checked.  At  "convergence",  7*  was  typically  within  a 
10~6  tolerance  of  7  (from  below).  Finally,  the  derivatives  with  respect  to  Ac,  Bc, 
and  Cc  were  checked  to  insure  they  were  close  to  zero  (they  should  be  exactly 
zero  for  the  true  solution).  Usually,  at  convergence,  most  of  the  elements  in 
these  derivative  matrices  were  around  10'3  to  10'5. 

5.2  Determining  an  Initial  Guess  for  DFP 

As  with  any  numerical  technique,  the  DFP  algorithm  requires  an  initial  guess 
of  the  solution.  As  discussed  in  Chapter  3,  the  central  Hw  (minimum  entropy) 
controller  is  the  solution  of  the  suboptimal  mixed  problem  for  n  =  1.  This  is  a 
full  order  compensator.  Therefore,  if  the  desired  order  of  the  mixed  solution  is 
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the  order  of  the  plant,  this  is  an  acceptable  initial  guess  for  the  given  level  of  7 
or  higher  (note  that  if  a  central  Hw  compensator  designed  for  a  certain  level  of 
7  is  used  as  a  starting  point  for  a  higher  7,  it  may  be  too  far  away  from  the 
solution  for  the  numerics  to  converge  even  though  it  is  an  "acceptable"  starting 
point).  While  the  full  order  minimum  entropy  controller  is  a  good  initial  guess 
for  a  full  order  mixed  problem,  it  cannot  be  used  directly  for  an  increased  order 
problem  because  its  size  is  incompatible.  Therefore,  finding  a  good  initial  guess 
for  the  increased  order  problem  is  not  quite  as  straightforward  as  the  full  order 
case. 

The  method  for  determining  an  initial  higher  order  guess  in  this  research 
utilized  the  J-Q  parameterization  of  K(s),  depicted  in  Figure  5-2. 


Figure  5-2.  Parameterization  of  K(s) 


If  K(s)  is  an  suboptimal  compensator,  J  is  completely  known.  The  central 
controller  (Q  =  0)  is  not  the  only  acceptable  choice  for  a  starting  guess.  As  long 
as  Q  is  stable,  strictly  proper,  and  has  ||  Q  ||  „  <  7,  it  will  produce  a 
compensator  that  admits  a  solution  to  the  X,  Y,  Q2,  and  Q*  Lyapunov  and 
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Riccati  equations  in  the  necessary  conditions.  Therefore,  the  initial  method  used 
was  to  simply  choose  a  Q  with  enough  arbitrary  left-hand  plane  poles  (and  an 
appropriate  gain  to  ensure  the  »-norm  bound  was  met)  to  give  the  desired  order 
of  K.  This  worked  for  the  n  +  1  order  compensator  in  the  SISO  example,  but  did 
not  work  for  any  higher  orders.  Even  though  the  compensator  was  theoretically 
acceptable,  it  was  not  close  enough  to  the  solution  for  the  algorithm  to  "start 
moving".  It  was  found  that  the  algorithm  is  very  sensitive  to  the  starting  guess. 

The  method  that  ended  up  being  the  most  reliable  was  less  arbitrary  than  the 
simplistic  approach  just  discussed.  It  is  possible  to  take  a  given  compensator  K 
and  a  J  (which  is  completely  specified  for  a  given  level  of  7)  to  calculate  a  Q 
that,  when  placed  in  a  feedback  loop  around  J,  will  produce  the  given  K  (see 
[Rid91a,152-154],  A  PRO-MATLAB™  routine  (named  K2Q )  that  accomplishes 
this  was  utilized.  It  was  found  that  the  Q  that  was  calculated  always  had  a 
relatively  high  order  (it  was  almost  never  a  minimal  realization). 

The  first  step  in  determining  a  higher  order  initial  guess  was  to  find  the 
minimum  entropy  controller  for  the  given  7.  If  this  compensator  is  input  into 
K2Q,  the  resulting  Q  would  be  Q  =  0  (as  it  should  be).  However,  if  the  central 
Hqo  compensator  was  modified  slightly  (by  simply  truncating  the  elements  in  its 
state  space  C  matrix  after  the  second  or  third  decimal  place)  the  resulting  Q 
would  be  non-zero  and  relatively  high  order.  Then,  using  a  balanced  Schur 
model  order  reduction,  Q  was  reduced  to  the  amount  of  states  (nc  -  n)  required 
to  achieve  the  desired  order  of  K.  Thus,  the  resulting  K  was  essentially  a 
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nonminimal  realization  of  the  central  H*  compensator.  This  process  does  not 
guarantee  that  the  ||  Q  ||  «  -  7  requirement  will  be  met.  Therefore,  the  state 
space  C  matrix  for  Q  would  be  arbitrarily  multiplied  by  0.6  in  order  to  reduce 
Q’s  maximum  singular  value.  This  was  always  done,  even  when  the  00-norm 
bound  on  Q  was  already  satisfied,  because  it  tended  to  help  the  algorithm  start 
moving. 

Once  a  mixed  solution  was  found,  it  could  also  be  used  as  a  starting  guess 
for  a  higher  compensator  order,  rather  than  the  minimum  entropy  controller.  Q 
was  always  non-zero  for  the  mixed  solution,  and  again  it  was  typically  high 
order.  The  same  process  of  reducing  Q  to  the  desired  order  and  multiplying  its 
C  matrix  by  0.6  was  also  used  for  coming  up  with  the  next  higher  order  starting 
point. 

Another  variable  that  has  not  been  discussed  yet  is  fi.  This  could  also  be 
varied  in  order  to  help  get  the  program  moving.  Usually,  the  first  initial  guess 
could  be  started  with  n  =  0.1;  however,  sometimes  n  had  to  be  increased  to  as 
much  as  0.7  in  order  to  start  the  process.  In  summary,  finding  a  good  initial 
guess  for  the  higher  order  case  is  not  a  trivial  task.  Many  times,  several 
different  attempts  had  to  be  made  before  an  acceptable  starting  point  could  be 
found. 
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5.3  SISO  Mixed  Optimization  Example 

Consider  the  mixed  H2/H00  optimization  system  block  diagram  that  was 
discussed  and  developed  in  Chapter  3  and  shown  in  Figure  5-3. 


Figure  5-3.  Mixed  Hj/H^  Optimization  Block  Diagram 


For  this  example,  all  signals  (d,  w,  e,  z,  u  and  y)  are  assumed  to  be  scalars. 
The  plant  P  has  the  state  space  realization 


r 

A  ! 

_ l 

Bd 

Bw 

Bu 

P(S)  - 

Ce  1 

i 

Ded 

B  ew 

Deu 

Cz 

Dzd 

B  ZW 

Dzu 

C  ' 

Dyd 

Byw 

D  yu 
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In  order  to  make  direct  comparisons  with  a  known  full  order  case,  the  same  SISO 
system  that  was  defined  in  [Rid9  la,  130-131]  was  chosen  for  this  analysis.  The 
state  space  matrices  of  the  system  are: 


A 


0.3908 

1.4453 

0.1288 


-0.4565 

-1.0491 

0.6744 


1.2657 

-1.2077 

1.0324 


0.0488 

1.4077 

-0.4275 

Bd  = 

0.3608 

0.3564 

K  = 

0.9723 

-1.6050 

Bu  - 

-0.4470 

-0.9172 

ce  =  [ 

0.9420 

0.0144 

cz  -  [ 

-0.0450 

0.3606 

cy  -  [ 

-1.5567 

-1.9432 

Ded 

=  [0] 

^ew  = 

Dzd 

«  [0] 

Dzw  = 

Dyd 

=  [0.5185] 

Dyw  = 

0.1187] 
1.8972  ] 
-0.0914  ] 


[0] 

^eu 

=  [  1.3575  ] 

[0] 

Dzu 

=  [  0.5781  ] 

[  0.3899  ] 

Dyu 

=  [0] 

First  of  all,  note  that  this  system  does  satisfy  all  the  assumptions  given  in  Section 
3.1  This  is  a  third  order  SISO  system  whose  "unweighted"  plant,  given  by 
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Pyu(s)  =  Cy(sI-A)'lBu  +  Dyu,  is  open-loop  unstable  and  minimum  phase.  The 
singular  value  plot  of  Pyu  is  shown  in  Figure  5-4. 


Figure  5-4.  Magnitude  Plot  of  SISO  Plant 


Before  beginning  the  mixed  problem,  it  will  be  helpful  to  know  the  limits  of 
achievable  H2  and  performance.  Consider  performing  pure  unconstrained  H2 
and  Hgo  optimization  on  the  given  plant.  Figure  5-5  shows  the  singular  value 
plot  of  the  unique  three-state  H2  optimal  controller  K2opt.  Figures  5-6  and  5-7 
are  the  singular  value  plots  of  the  corresponding  closed-loop  transfer  functions 
Tzw  and  Ted.  The  minimum  achievable  2-norm  of  Tzw  is  ac  =  9.9263.  The 
oo-norm  of  Ted  using  K2opt  is  y2  =  4.5364. 
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Figure  5-5.  Singular  Value  Plot  of  K2opt 


frequency  (rad/s) 


Figure  5-6.  Singular  Value  Plot  of  Tw  for  K2opt 
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Figure  5-7.  Singular  Value  Plot  of  Ted  for  K2opt 


Now,  if  optimization  is  performed,  the  minimum  achievable  oo-norm  of  Ted 

is  found  to  be  about  70  *  2.1426.  The  freedom  parameter  Q(s)  from  the 
parameterization  of  controllers  must  be  specified.  The  central 
compensator  (i.e.  minimum  entropy  controller)  is  of  particular  interest  since  it 
will  be  the  initial  guess  for  the  DFP  program  (at  n  =  1),  so  choose  Q(s)  =  0. 
The  singular  value  plot  of  the  Hw  suboptimal  central  compensator  for  7  =  2.1426 
(Koo2.i426)  ls  giyen  in  Figure  5-8.  Figures  5-9  and  5-10  are  the  singular  value 
plots  of  the  corresponding  closed-loop  transfer  functions  Tzw  and  Ted.  The 
2-norm  of  Tzw  for  K,^  l426  is  ||Tzw||2  =  94.323,  considerably  higher  than  a0. 
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frequency  (rad/s) 


Figure  5-8.  Singular  Value  Plot  of  £<,,,2.1426 


frequency  (rad/s) 

Figure  5-9.  Singular  Value  Plot  of  Tw  for  ^2.1426 
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T3 

a  2.15- 

a  _ 

£  2.14- 

2.13- 

2.12- 

2.11- 

2.  £  _ > _ 1  1  1  1  1 1 U  ■  i  i  i  i  i  in  i  i  i  i_j_i  m  i  i  i  i  i  i  in - i  i  i  i  i  i  in  i  i  ■  i  i  i  n 
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frequency  (rad/s) 

Figure  5-10.  Singular  Value  Plot  of  Ted  for  Koo2.i426 

Although  it  is  not  evident  in  these  plots,  they  all  have  a  high  frequency  roll  off. 
If  the  optimal  HM  controller  was  found,  there  would  be  no  roll  off  and  ||Tzw||2 
would  be  infinite. 

Now  consider  performing  mixed  Hi/H,*  optimization  assuming  the  controller 
is  full  order.  A  brief  summary  of  Ridgely’s  results  are  shown  in  Table  5-1  and 
Figure  5-11.  Table  5- 1  shows  the  values  of  ||  Tzw  ||  2  and  |j  Ted  ||  „  for  the  mixed 
and  central  H,*  controllers  at  varying  levels  of  y.  Notice  that  ||Ted||00=7  for 
the  mixed  controller.  Actually,  the  oo-norm  values  are  rounded  off  (typically 
within  0.00001).  Recall  that  the  true  mixed  solution  must  lie  on  the  boundary 
of  the  oo-norm  constraint.  Figure  5-11  shows  this  data  in  graphical  form. 
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Table  5-1.  SISO  Example  Full  Order  Results 


(Q  =  0) 

O 

II 

o 

(mix) 

(mix) 

7 

1  Te<1  II . 

II  Tzw  II  2 

l|TedIL 

l|T2J2 

2.1426 

2.1426 

94.323 

2.1426 

~  94.323 

2.145 

2.145 

20.387 

2.145 

~  20.387 

2.15 

2.15 

16.099 

2.15 

16.091 

2.2 

2.2 

14.055 

2.2 

13.868 

2.25 

2.247 

14.076 

2.25 

13.295 

2.35 

2.3389 

14.305 

2.35 

12.501 

2.5 

2.4675 

14.622 

2.5 

11.616 

2.75 

2.6589 

14.953 

2.75 

10.870 

3.0 

2.8243 

15.107 

3.0 

10.460 

3.25 

2.9673 

15.161 

3.25 

10.225 

3.5 

3.0909 

15.164 

3.5 

10.086 

4.0 

3.2913 

15.107 

4.0 

9.9539 

4.5364 

3.4532 

4.5364 

9.9263 

10.0 

3.9948 

14.666 

4.5364 

9.9263 

50.0 

4.1559 

14.557 

4.5364 

9.9263 

100.0 

4.1611 

14.554 

4.5364 

9.9263 

[Rid91a,  166] 
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Figure  5-11.  SISO  Example  Full  Order  and  Q(s)=0  Results 

[Rid91a,167] 

No  discussion  of  these  full  order  results  will  be  given  here.  Rather,  these  results 
are  included  as  a  point  of  reference  for  the  higher  order  results. 

5.3.1  7  =  2.5  Results.  Now,  with  all  the  preliminaries  taken  care  of, 

consider  the  case  of  higher  order  compensators.  Since  the  nature  of  the  results 
using  a  higher  order  compensator  was  unknown,  the  basic  approach  was  to  begin 
with  the  full  order  case  and  continually  increase  the  compensator  order  until  some 
kind  of  trend  could  be  recognized.  Before  this  order  sweep  could  be 
accomplished,  the  design  7  had  to  be  chosen.  7  =  2.5  was  the  first  level  selected. 
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Then,  compensators  of  the  following  orders  were  obtained  by  running  DFP:  3, 
4,  5,  6,  7,  8,  9,  12,  18.  Table  5-2  shows  a  summary  of  the  higher  order  results 
for  7  =  2.5.  Note  that  the  »-norm  values  are  essentially  the  same  as  y.  They 
are  actually  slightly  less  than  y  (typically  by  about  10'6). 


Table  5-2.  Higher  Order  Results,  7  =  2.5 


Compensator 

Order 

* 

a 

* 

7 

3  (full) 

11.61625 

2.5 

4 

11.56329 

2.5 

5 

11.52972 

2.5 

6 

11.51510 

2.5 

7 

11.50500 

2.5 

8 

11.49373 

2.5 

9 

11.48754 

2.5 

12 

11.48654 

2.5 

18 

11.48647 

2.5 

Figure  5-12  shows  this  same  data  on  a  graph.  Note  that  since  only  integer  order 
compensators  are  allowed,  this  is  not  a  continuous  curve.  In  [Rid91a]  it  was 
shown  that  a*  =  1 1.507  for  nc  =  9  at  7  =  2.5  (this  was  the  only  higher  order  result 
given).  The  slight  discrepancy  is  easily  accounted  for  by  the  fact  that  the 
question  of  "close  enough"  to  the  true  solution  is  fairly  arbitrary.  Apparently, 
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the  solution  given  here  is  converged  a  little  more.  In  fact,  it  is  possible  to  very 


slightly  reduce  the  values  given  here  even  more. 


11.62 


11.6 


11.58 

t  11-56 

o 
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°  11.54 

I 

CM 
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11.5 

11.48 

2  4  6  8  10  12  14  16  18  20 

Compensator  Order 

Figure  5-12.  Higher  Order  Results,  7=2.5 
One  observation  can  immediately  be  made:  in  general,  the  optimal  order  is 


definitely  not  the  order  of  the  plant.  This  graph  clearly  shows  that  as  order  is 


increased  beyond  full  order  (n  =  3),  the  value  of  the  2-norm  of  Tzw  decreases. 


There  does  appear  to  be  an  order  beyond  which  meaningful  reductions  in  a*  are 


no  longer  attained  (i.e.  nc  =  9).  Also,  in  this  example,  a*  is  strictly 


monotonically  decreasing  with  increasing  order.  This  behavior  has  not  been 


proven  in  general  (if  it  was,  it  would  prove  that  the  optimal  order  is  infinite). 


However,  it  does  appear  to  be  true  in  this  example. 
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Figures  5-13,  5-14,  and  5-15  show  the  singular  value  plots  of  the  mixed  ItyH,* 
compensators  with  increasingly  higher  orders.  Notice  first  that  they  all  have  the 
same  basic  shape.  The  sharp  peak  at  about  2.5  rad/sec  decreases  in  magnitude 
with  increasing  order.  Also,  as  might  be  expected  from  Figure  5-12,  there  does 
not  appear  to  be  much  difference  between  the  9,  12,  and  18-state  compensators 
(at  least  in  terms  of  frequency  response).  Whether  or  not  these  compensators  are 
actually  different  realizations  of  the  same  9th  order  controller  will  be  discussed 
later.  From  these  plots,  they  certainly  appear  to  be  essentially  the  same. 


frequency  (rad/s) 

Figure  5-13.  Singular  Value  Plots  of  (3,4,5,6-state) 
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magnitude  magnitude 


frequency  (rad/s) 

Figure  5-14.  Singular  Value  Plots  of  (6,7,8,9-state) 


frequency  (rad/s) 

Figure  5-15.  Singular  Value  Plots  of  (9,12,18-state) 
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Figures  5-16,  5-17,  and  5-18  show  the  singular  value  plots  of  Tzw  for  the  higher 
order  mixed  solutions.  In  these  plots  it  is  apparent  that  the  higher  order 
compensators  decrease  the  2-norm  of  Tzw  by  making  the  peak  in  the  curve 
progressively  narrower.  Even  though  the  peak  rises  slightly,  the  net  change  in 
area  under  the  curve  decreases.  Again,  there  are  no  distinguishable  differences 
in  the  9,  12  and  18-state  results.  Figures  5-19,  5-20,  and  5-21  show  the  singular 
value  plots  of  Ted  for  the  higher  order  mixed  solutions.  These  plots  are  the  most 
interesting  because  they  best  demonstrate  the  value  of  the  higher  order 
controllers.  In  [Rid9 la] ,  it  was  shown  for  the  full  order  case  that  the  mixed 
solution  tries  to  recover  to  the  H2  optimal  solution.  This  makes  sense  since  both 
problems  have  the  same  performance  objective.  However,  since  the  mixed 
problem  has  a  constraint  that  must  be  met,  it  will  try  to  match  the  true  optimal 
solution  as  best  as  it  can  while  satisfying  the  constraint.  This  can  be  seen 
dramatically  in  Figure  5-22.  This  plot  shows  the  H2  optimal,  the  central  H^,, 
and  the  full  order  mixed  solutions.  Notice  how  the  mixed  solution  tries  to 
recover  to  the  H2  solution.  In  the  regions  where  the  H2  curve  exceeds  the 
oo -norm  bound,  the  mixed  solution  basically  lies  right  on  the  oo-norm  boundary. 
Figure  5-23  shows  an  expanded  view  of  this  same  plot  with  the  9-state  solution 
included  (and  Q=0  omitted).  The  extra  degrees  of  freedom  provided  by  the 
higher  order  compensator  enable  a  better  recovery  of  the  H2  solution.  As  can  be 
seen,  the  9-state  curve  makes  a  much  sharper  turn  where  the  H2  and  Q  =  0  curves 
intersect.  It  also  dips  further  down  into  the  notch. 


magnitude  magnitude 


10-3  10-2  10-1  10°  10‘  102  103 

frequency  (rad/s) 

Figure  5-16.  Singular  Value  Plots  of  Tw  for  7=2.5 
(3,4,5,6-state) 
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Figure  5-17.  Singular  Value  Plots  of  Tw  for  7=2.5 
(6,7, 8, 9- state) 
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magnitude  magnitude 


Figure  5-18.  Singular  Value  Plots  of  Tw  for  7=2.5 
(9,12,18-state) 


frequency  (rad/s) 

Figure  5-19.  Singular  Value  Plots  of  Ted  for  7=2.5 
(3,4,5,6-state) 
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magnitude 


frequency  (rad/s) 

Figure  5-20.  Singular  Value  Plots  of  Ted  for  7=2.5 
(6,7,8,9-state) 


Figure  5-21.  Singular  Value  Plots  of  Ted  for  7=2.5 
(9,12,18-state) 
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l0-3  10-a  10-1  10°  10*  10*  103 

frequency  (rad/s) 

Figure  5-22.  Ted  Comparison  Plot  for  7=2.5 


frequency  (rad/s) 

Figure  5-23.  Ted  Comparison  Plot  Expanded  View  (7=2.5) 


5-26 


From  all  these  plots,  it  is  clear  that  increasing  compensator  order  does  produce 
a  better  solution.  However,  after  9th  order,  no  more  significant  improvement  is 
seen.  The  question  then  becomes:  is  3n  an  optimal  order?  Are  the  12th  and  18th 
order  controllers  simply  nonminimal  realizations  of  the  9th  order  controller? 
Several  different  approaches  were  taken  to  answer  this  question.  The  first 
approach  was  to  examine  the  Hankel  singular  values  of  the  compensators.  If  the 
9-state  controller  truly  has  an  optimal  order,  it  might  be  expected  that  the  9-state 
controller  would  have  9  relatively  signficant  Hankel  singular  values.  The 
12-state  controller  would  have  3  Hankel  singular  values  that  are  either  exactly 
zero  or  at  least  several  orders  of  magnitude  smaller  than  the  other  9.  Likewise, 
the  18-state  contoller  would  have  9  relatively  small  Hankel  singular  values. 
However,  this  phenomenon  was  not  found  to  be  true.  Table  5-3  shows  a 
summary  of  the  Hankel  singular  values  of  Kraix  for  each  order,  in  rough 
magnitudes.  Notice  that  even  though  some  of  the  Hankel  singular  values  get 
relatively  small  for  the  higher  order  compensators,  there  is  never  any  definite 
division  where  the  remainder  of  the  states  are  clearly  superfluous.  Thus,  at  first 
glance,  this  approach  does  not  give  much  help  in  determining  if  the  higher  order 
compensators  should  be  reduced. 

Next,  the  actual  poles  and  zeros  of  the  compensators  were  plotted  in  order 
to  get  a  different  look  at  possible  pole/zero  cancellations.  Unfortunately,  this  did 
not  provide  any  more  help.  Every  compensator  6th  order  and  up  had  multiple 
poles  with  zeros  right  on  top  of  them.  Simply  by  inspection,  if  the  approximate 
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Table  5-3.  Hankel  Singular  Values  of  Kmix  (7=2.5) 


3 

state 

4 

state 

5 

state 

6 

state 

7 

state 

8 

state 

9 

state 

12 

state 

18 

state 

7e-0 

6e-0 

2e-0 

7e-0 

6e-0 

2e-0 

le-1 

7e-0 

6e-0 

2e-0 

2e-l 

3e-2 

6e-0 

6e-0 

2e-0 

le-1 

4e-2 

3e-2 

7e-0 

6e-0 

2e-0 

le-1 

7e-2 

4e-2 

3e-3 

3e-0 

3e-0 

3e-0 

3e-0 

2e-0 

4e-l 

5e-2 

4e-2 

4e-0 

4e-0 

3e-0 

3e-0 

2e-0 

le-1 

5e-2 

4e-2 

2e-3 

6e-0 

5e-0 

3e-0 

le-0 

9e-l 

2e-l 

7e-2 

5e-2 

4e-2 

2e-2 

9e-3 

3e-4 

5e-0 

5e-0 

5e-0 

3e-0 

2e-0 

le-1 

5e-2 

4e-2 

le-2 

le-2 

6e-3 

4e-4 

2e-4 

7e-5 

3e-5 

3e-5 

4e-6 

4e-7 

pole/zero  cancellations  are  made,  each  of  the  higher  order  compensators  would 
reduce  to 


6-state 

4-state 

7- state 

4-state 

8-state 

—¥ 

4 -state 

9-state 

- 

(4  or  6)-state 

12-state 

- 

(7  or  8)-state 

18-state 

- 

(8  or  10)-state 
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The  final  approach  that  was  tried  was  to  simply  perform  a  Schur  balanced 
model  order  reduction  on  the  12  and  18-state  compensators  to  see  what  they 
looked  like  when  reduced  to  9th  order.  It  was  quickly  discovered,  however,  that 
it  is  impossible  to  remove  even  one  state  without  immediately  violating  the 
oo -norm  bound  on  Ted.  This  means  that  these  reduced  order  controllers  are  no 
longer  valid  solutions  to  the  mixed  problem.  It  is  not  really  surprising  that  this 
happens  because  the  solution  (for  each  order)  lies  right  on  the  boundary  of  the 
oo -norm  constraint.  Therefore,  any  reduction  in  order  immediately  makes  the 
solution  unacceptable.  Based  on  this  result,  it  seems  that  all  the  apparent 
pole/zero  cancellations  are  not  true  cancellations.  These  additional  poles  may 
have  very  small  residues,  but  they  need  to  remain  in  the  solution. 


5.3.2  7  =  3.0  Results.  In  order  to  see  if  the  same  trends  observed  in  the 
7=2.5  case  are  typical,  the  7  level  was  fixed  at  3.0.  This  time,  however,  only 
6,  9,  and  12-state  solutions  were  found.  Table  5-4  shows  a  summary  of  the 
higher  order  results  for  7  =  3.0. 


Table  5-4.  Higher  Order  Results,  7  =  3.0 


Compensator 

Order 

* 

a 

* 

7 

3  (full) 

10.460 

3.0 

6 

10.44380 

3.0 

9 

10.42787 

3.0 

12 

10.42619 

3.0 

Figure  5-24  shows  this  same  data  graphically. 
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Figure  5-24.  Higher  Order  Results,  7=3.0 

The  same  decrease  in  a*  for  increased  order  compensators  is  clearly  seen  for  this 
level  of  7  also.  Notice,  however,  that  there  is  less  of  an  improvement  than  for 
the  7=2.5  case.  This  is  because  the  full  order  mixed  solution  is  already  closer 
to  the  H2  optimal  solution  due  to  the  relaxation  of  7.  Figures  5-25,  5-26,  and 
5-27  show  the  singular  value  plots  of  Kmix,  Tzw,  and  Ted  for  orders  6,  9  and  12. 
The  same  trends  that  were  observed  for  7  =  2.5  are  apparent  here.  The  recovery 
of  the  H2  solution  is  again  seen.  As  shown  in  Figure  5-28,  the  7  =  3.0  solution 
is  able  to  recover  more  of  the  H2  solution  due  to  the  higher  level  of  7.  Also,  the 
higher  order  compensators  do  a  better  job  with  this  recovery. 
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magnitude  magnitude 


Figure  5-25.  Singular  Value  Plot  of  K,^,  7=3.0 
(6,9,12-state) 


Figure  5-26.  Singular  Value  Plot  of  T™,  7=3.0 
(6,9,12-state) 

5-31 


magnitude 


Figure  5-27.  Singular  Value  Plot  of  Ted,  7=3.0 
(6,9, 12-state) 


frequency  (rad/s) 

Figure  5-28.  Ted  Comparison  Plot  for  7=3.0 
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5.3.3  Ninth-Order  Results,  y  Sweep.  Now,  consider  freezing  the 
compensator  order  and  finding  the  solutions  for  a  whole  range  of  7’s.  Since  the 
9-state  controllers  for  the  two  7  levels  examined  so  far  appear  to  be  near  the 
limit  of  achievable  mixed  performance,  select  the  compensator  order  to  be  9. 
Table  5-5  shows  a  summary  of  the  results. 


Figure  5-5.  Results  of  9-state  7  Sweep 


7 

3-state  Kmix 
a 

9-state  Kmix 
a 

Percent 

Difference 

2.2 

13.868 

13.697 

1.230 

2.25 

13.295 

13.085 

1.578 

2.35 

12.501 

12.243 

2.060 

2.5 

11.616 

11.488 

1.105 

2.75 

10.870 

10.795 

0.690 

3.0 

10.460 

10.428 

0.307 

3.25 

10.225 

10.219 

0.059 

3.5 

10.086 

10.078 

0.079 

Note  that  the  3-state  solutions  are  taken  directly  from  Table  5-1.  Figure  5-29 
shows  Ridgley’s  full  order  mixed  plot  with  the  9th  order  solutions  superimposed. 
Figure  5-30  is  an  expanded  view  of  this  same  plot.  One  fact  is  very  obvious: 
the  9th  order  compensators  do  not  give  large  improvements  in  terms  of  reducing 
the  2-norm  for  this  example.  The  largest  decrease  in  a*  is  about  2%. 
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-norm  of  Tzw,  mix  2-norm  of  Tzw,  mix 


20 


18 


Inf-norm  of  Ted,  mix 

Figure  5-30.  Mixed  Plot  (expanded),  3  and  9-state  Comparison 


5-34 


The  mixed  solution  singular  value  plots  of  Ted  and  Tzw  for  the  9th  order 
compensators  are  given  in  Figures  5-31  and  5-32.  These  plots  show  a  definite 
recovery  of  the  H2  optimal  solution.  This  was  true  for  the  full  order  case 
[Rid91a,164],  and  is  seen  to  be  true  for  the  higher  order  case  as  well.  In  fact, 
as  shown  earlier,  the  higher  order  compensators  actually  do  a  better  job  of 
recovering  the  unconstrained  H2  optimal  solution  than  the  full  order  controllers. 


frequency  (rad/s) 

Figure  5-31.  Mixed  Ted  Plot  (9-state) 

7  =  2.1426,  2.25,  2.5,  2.75,  3.0,  3.5,  4.5364 


frequency  (rad/s) 

Figure  5-32.  Mixed  Tw  Plot  (9-state) 

7  =  2.1426,  2.25,  2.5,  2.75,  3.0,  3.5,  4.5364 
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5.4  MIMO  Mixed  Optimization  Example 


Consider  the  mixed  H2/H0o  optimization  system  block  diagram  shown  in 
Figure  5-3.  For  this  multi-input  multi-output  example,  all  signals  (d,  w,  e,  z, 
u  and  y)  are  assumed  to  be  two-dimensional  vectors.  In  order  to  make  direct 
comparisons  with  a  known  full  order  case,  the  same  system  that  was  defined  in 
[Rid9 la,  170-171]  was  chosen  for  this  analysis.  The  state  space  matrices  of  the 
system  are: 


-5  2  14  20 

10  0  0 
0  10  0 

0  0  10 


0.03 

0.008 

0.22 

0.93 

0.07 

0.44 

0.05 

0.38 

= 

0.05 

0.38 

B„  = 

0.63 

0.77 

0.53 

0.07 

W 

0.68 

0.52 

U 

0.88 

0.48 

0.67 

0.42 

0.58 

0.83 

0.27 

0.24 

0.55  0.33  1.80  0.12 

0.72  0.97  1.82  1.81 


0.07  0.38  0.91  0.46 
0.50  0.28  0.53  0.94 

0.05  0.77  0.13  0.69 
0.76  0.83  0.02  0.87 
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Note  that  this  system  does  satisfy  all  the  assumptions  given  in  Section  3.1  This 
is  a  fourth  order  MIMO  system  whose  "unweighted"  plant,  given  by 
Pyu(s)  =  Cyfsl-A)'1!^  +  Dyu,  is  open-loop  unstable  and  nonminimum  phase. 
The  singular  value  plot  of  P  is  shown  in  Figure  5-33. 


frequency  (rad/s) 

Figure  5-33.  Singular  Value  Plot  of  MIMO  Plant 


Before  beginning  the  mixed  problem,  it  will  be  helpful  to  know  the  limits  of 
achievable  H2  and  Ha  performance.  Consider  performing  pure  unconstrained  H2 
and  Hm  optimization  on  the  given  plant.  Figure  5-34  shows  the  singular  value 
plot  of  the  unique  four-state  H2  optimal  controller  K2opt.  Figures  5-35  and  5-36 
are  the  singular  value  plots  of  the  corresponding  closed-loop  transfer  functions 
Tzw  and  Ted-  The  minimum  achievable  2-norm  of  Tzw  is  aQ  =  0.7975.  The 
oo-norm  of  Ted  using  K2opt  is  y2  =  40.548. 


frequency  (rad/s) 

Figure  5-34.  Singular  Value  Plot  of  K2opt 
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magnitude  magnitude 


5 


frequency  (rad/s) 


Figure  5-35.  Singular  Value  Plot  of  Tzw  for  K2opt 


frequency  (rad/s) 

Figure  5-36.  Singular  Value  Plot  of  Ted  for  K2opt 
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Now,  if  H^,  optimization  is  performed,  the  minimum  achievable  oo-norm  of  Ted 
is  found  to  be  about  70  *  2.3012.  The  freedom  parameter  Q(s)  from  the 
parameterization  of  H^,  controllers  must  be  specified.  The  central  Ha 
compensator  is  of  particular  interest  since  it  will  be  the  initial  guess  for  the  DFP 
program  (at  n  =  1),  so  choose  Q(s)  =  0.  The  singular  value  plot  of  the 
suboptimal  central  compensator  for  7=2.3012  (K^^on)  is  given  in  Figure  5-37. 
Figures  5-38  and  5-39  are  the  singular  value  plots  of  the  corresponding  closed- 
loop  transfer  functions  Tzw  and  Ted.  The  2-norm  of  Tzw  for  1^2.3012  is 
| Tzw  1 2  =  1887.3,  considerably  higher  than  aQ. 


frequency  (rad/s) 

Figure  5-37.  Singular  Value  Plot  of  Koo2.3oi2 
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magnitude  magnitude 


frequency  (rad/s) 


Figure  5-38.  Singular  Value  Plot  of  Tw  for  1^2.3012 


frequency  (rad/s) 

Figure  5-39.  Singular  Value  Plot  of  Ted  for  1^2.3012 


5-42 


Although  it  is  not  evident  in  these  plots,  they  all  have  a  high  frequency  roll  off. 
If  the  optimal  controller  was  found,  there  would  be  no  roll  off  and  |Tzw|2 
would  be  infinite. 

Now  consider  performing  mixed  optimization,  assuming  the  controller 

is  full  order.  A  brief  summary  of  some  of  Ridgely’s  MIMO  results  are  shown 
in  Table  5-6  and  Figure  5-40.  Table  5-6  shows  the  values  of  ||  Tzw  ||  2  and 
II  Te(j  II  oo  f°r  the  mixed  and  central  HM  controllers  at  varying  levels  of  7.  Notice 
that  1  Ted  ||  oo  =  7  for  the  mixed  controller.  Actually,  the  00 -norm  values  are 
rounded  off  (typically  within  0.00001).  Figure  5-40  shows  this  data  in  graphical 
form.  Note  that  y2  is  not  shown  on  the  plot  because  it  is  out  at  40.548.  This 
shows  the  tremendous  value  of  the  mixed  controllers.  The  long,  flat  "tail"  on  the 
mixed  curve  enables  huge  reductions  in  the  00 -norm  of  Ted  while  making  only 
small  sacrifices  in  the  2-norm  of  Tzw.  No  further  discussion  of  these  full  order 
results  will  be  given  here.  Rather,  these  results  are  included  as  a  point  of 
reference  for  the  higher  order  results. 
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Table  5-6.  MIMO  Example  Full  Order  Results 


(Q  =  0) 

— \ 
o 

II 

a 

(mix) 

(mix) 

y 

IITJL 

!ITzw»2 

IITJL 

|TIW||2 

2.3012 

2.3012 

1887.3 

2.3012 

~  1887.3 

2.32 

2.3199 

119.129 

2.32 

115.713 

2.35 

2.3495 

62.307 

2.35 

47.619 

2.4 

2.3979 

37.008 

2.4 

23.664 

2.5 

2.4914 

22.155 

2.5 

10.713 

2.6 

2.5807 

16.799 

2.6 

6.9406 

2.75 

2.7069 

13.11 

2.75 

4.6072 

3.0 

2.8975 

10.415 

3.0 

3.1622 

4.0 

3.4578 

7.5491 

4.0 

1.9115 

5.0 

3.7978 

6.8949 

5.0 

1.4582 

10.0 

4.3708 

6.2236 

10.0 

0.9800 

40.548 

4.5874 

6.0753 

40.548 

0.7975 

50.0 

4.5925 

6.0723 

40.548 

0.7975 

100.0 

4.5997 

6.0678 

40.548 

0.7975 

[Rid9 la, 1 88] 
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Figure  5-40.  MIMO  Example  Full  Order  and  Q(s)=0  Results 


[Rid91a,189] 

Now,  consider  the  case  of  higher  order  compensators.  As  in  the  SISO 
example,  the  basic  approach  in  the  MIMO  problem  was  to  begin  with  the  full 
order  case  and  continually  increase  the  compensator  order  until  some  kind  of 
trend  could  be  recognized.  Before  this  order  sweep  could  be  accomplished,  the 
design  7  had  to  be  chosen.  As  7  approaches  70,  the  numerics  become  more  and 
more  difficult.  Therefore,  a  level  of  7  =  3.0  was  selected.  This  was  the  only  7 
that  was  run  due  to  the  excessive  computer  time  required  to  find  a  solution  (for 
example,  the  entries  in  Table  5-7  typically  took  from  about  eight  to  twelve  hours 
for  the  4-state  solution  to  about  two  weeks  of  almost  continuous  run-time  for  the 
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16-state  solution).  In  retrospect,  perhaps  a  slightly  lower  y  might  have  been  a 
better  choice  to  demonstrate  the  value  of  the  higher  order  compensators.  As  it 
turned  out,  the  improvements  due  to  the  higher  order  controllers  were 
observable,  but  they  were  small.  Compensators  of  the  following  orders  were 
then  obtained  by  running  DFP:  4,  6,  7,  8,  9,  10,  12,  16. 

Table  5-7  shows  a  summary  of  the  higher  order  results  for  y  =  2.5.  Note 
that  the  oo-norm  values  are  essentially  the  same  as  y.  They  are  actually  slightly 
less  than  y  (typically  by  about  10'6). 


Table  5-7.  Higher  Order  Results,  7  =  3.0 


Compensator 

Order 

* 

a 

* 

7 

4  (full) 

3.15799 

3.0 

6 

3.15530 

3.0 

7 

3.15411 

3.0 

8 

3.15246 

3.0 

9 

3.15238 

3.0 

10 

3.15232 

3.0 

12 

3.15231 

3.0 

16 

3.15216 

3.0 

Figure  5-12  shows  this  same  data  on  a  graph. 
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Figure  5-41.  Higher  Order  Results,  7=3.0 

One  observation  can  immediately  be  made:  in  general,  the  optimal  order  is 
definitely  not  the  order  of  the  plant.  This  graph  clearly  shows  that  as  order  is 
increased  beyond  full  order  (n=4),  the  value  of  the  2-norm  of  Tzw  decreases. 
This  plot  does  not  have  a  nice  exponential-type  decay  as  seen  in  the  SISO 
example.  It  might  be  that  some  of  the  solutions  are  not  completely  converged. 
On  the  other  hand,  there  is  nothing  that  says  what  the  shape  of  this  curve  must 
be.  In  fact,  even  though  this  curve  is  also  strictly  monotonically  decreasing, 
there  are  no  known  proofs  that  guarantee  this  in  general.  The  key  observation 
here  is  that  all  the  higher  order  compensators  do  have  a  lower  a*  than  the  full 
order  solution.  Also,  notice  that  this  curve  basically  bottoms  out  at  nc  =  8. 
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This  is  two  times  the  order  of  the  plant  Recall  that  in  the  SISO  example,  this 
happened  at  three  times  the  order  of  the  plant). 

The  singular  value  plots  of  the  mixed  solutions  are  given  in  Figures  5-42  and 
5-43.  The  corresponding  singular  value  plots  of  Tzw  are  shown  in  Figures  5-44 
and  5-45.  Figures  5-46  and  5-47  are  the  corresponding  singular  value  plots  of 
Ted.  The  same  trends  that  were  observed  in  the  SISO  case  can  be  seen  here  in 
the  MIMO  case.  As  already  mentioned,  the  higher  order  compensators  do  not 
produce  large  changes  or  improvements  in  this  example.  This  is  probably  a 
combination  of  the  system  itself  and  the  choice  of  7.  However,  even  though  the 
changes  are  small,  they  do  produce  some  improvement.  In  order  to  show  the 
recovery  of  the  mixed  solution  to  the  H2  solution,  Figure  5-48  shows  an 
expanded  view  of  the  Ted  plots.  The  H2  solution,  4-state  mixed,  and  16-state 
mixed  solutions  are  given  for  comparison.  Notice  that  even  though  the  mixed 
solutions  do  not  make  the  sharp  turn  to  follow  the  H2  curve  immediately,  they 
do  make  the  sharp  turn.  The  16-state  mixed  compensator  makes  this  turn  sharper 
than  the  4-state  compensator  (as  seen  in  the  SISO  example). 
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magnitude  magnitude 


frequency  (rad/s) 

Figure  5-43.  Singular  Value  Plots  of  (9,10,12,16-state) 
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ma«nitude  magnitude 


Figure  5-44.  Singular  Value  Plots  of  Tw  for  7=3.0 
(4,6,7,8-state) 


Figure  5-45.  Singular  Value  Plots  of  Tw  for  7=3.0 
(9,10,12,16-state) 
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Figure  5-46.  Singular  Value  Plots  of  Ted  for  7=3.0 
(4,6,7,8-state) 


frequency  (rad/s) 

Figure  5-47.  Singular  Value  Plots  of  Ted  for  7=3.0 
(9,10,12,16-state) 
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VI.  Conclusions  and  Recommendations 

6.1  Optimal  Order 

Extensions  from  the  full  order  case  to  the  increased  order  case  were  made 
in  Chapter  IV,  and  it  has  been  shown  by  numerous  examples  that  the  optimal 
order  of  the  mixed  solution  is  not,  in  general,  the  order  of  the  plant. 

This  is  a  departure  from  the  nature  of  the  separate,  unconstrained  H2  and  Ha 
problems,  and  only  serves  to  demonstrate  the  complexity  of  the  mixed  problem. 
While  it  is  now  known  that  the  optimal  order  of  the  mixed  problem  is  not,  in 
general,  the  order  of  the  plant,  the  ultimate  question  of  what  the  optimal  order 
is  remains  to  be  proven. 

At  the  outset,  it  was  believed  that  the  optimal  order  might  be  3n  (and  some 
of  the  initial  findings  seemed  to  confirm  this.)  A  brief  outline  of  the  rationale 
for  this  conjecture  is  as  follows.  Since  the  oo-norm  bound  must  be  satisfied,  the 
compensator  K(s)  can  be  parameterized  as  a  lower  LFT  of  J  and  Q  as  shown  in 
Figure  6-1,  where  J  is  given  by  an  HB  parameterization  [DGKF89]  and 
Q€RH„,  |  Q  |  *  <7. 
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Figure  6-1,  Mixed  Optimization  Block  Diagram  with  Q(s) 


Since  J  is  completely  known,  it  can  be  combined  with  the  nominal  plant  P  to 
form  a  new  plant  P}  whose  order  is  2n  (assuming  no  pole-zero  cancellations 
occur).  This  is  shown  in  Figure  6-2. 


Figure  6-2.  Closure  of  the  P-J  Loop  Through  Optimization 


The  problem  is  now 

inf  II  Tzw  II 2  subject  to  1 Q  |  a,  <  7 
QERHge 

If  it  can  be  shown  that  the  optimal  order  of  Q  is  the  order  of  this  new  plant  (2n), 
the  resulting  compensator  would  have  order  3n.  Unfortunately,  the  problem  set 
up  in  this  manner  is  no  more  tractable  than  the  original.  Therefore,  no  analytical 
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solution  to  this  problem  has  been  found  either  (if  it  had,  the  general  mixed 
problem  would  be  solved).  Also,  this  Q  could  be  further  parameterized  by  a  J2 
and  Q2  from  an  H2  parameterization  on  Pj,  and  a  new  "plant"  of  order  3n 
created.  If  it  could  be  shown  that  this  continues  indefinitely,  it  would  lead  to  an 
optimal  mixed  compensator  of  infinite  order. 

Since  it  has  been  shown  that  adding  states  to  the  controller  does  cause 
improvement,  it  seems  intuitive  that  this  may  continue  indefinitely  while 
asymptotically  approaching  some  minimum  achievable  a*.  In  fact,  as  a 
conclusion  of  this  research,  this  is  formally  conjectured  and  supported  with  three 
evidences. 

Conjecture:  Assume  7  is  selected  such  that  yQ  <  7  <  72.  The  optimal  order 
of  Kmix  is  infinite. 

Evidence  1:  In  both  the  SISO  and  MIMO  examples,  every  increase  in  order  of 
Kraix  produced  a  reduction  in  the  2-norm  of  Tzw.  Now,  it  could  be  rightfully 
argued  that  these  reductions  may  be  within  the  noise  level  of  the  computational 
abilities  of  the  computer.  Therefore,  the  reductions  in  the  2-norm  are  not  real. 
However,  care  was  taken  to  perform  all  calculations  in  double  precision.  Also, 
it  does  seem  significant  that  reductions  occurred  in  every  instance.  The  numbers 
being  dealt  with  may  be  small,  but  they  should  not  be  completely  discounted. 
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Evidence  2:  It  has  been  demonstrated  that  the  mixed  solution  tries  to  recover  to 


the  H2  solution  while  still  meeting  the  oo-norm  constraint.  This  is  seen  most 
clearly  on  the  Ted  singular  value  plots.  Figure  5-23  is  repeated  below  in  Figure 
6-3  for  the  sake  of  discussion. 


frequency  (rad/s) 

Figure  6-3.  SISO  Ted  Comparison  Plot  Expanded  View  (7=2.5) 


The  mixed  solution  lies  on  the  oo-norm  boundary  wherever  the  H2  solution  is 
above  the  design  y  level.  When  the  H2  solution  drops  below  the  oo-norm  bound, 
the  mixed  solution  tries  to  follow  this  curve.  Upon  inspection  of  the  higher  order 
results,  it  appears  that  the  point  where  the  mixed  curve  turns  to  follow  the  H2 
curve,  the  true  optimal  mixed  solution  may  have  a  point  of  discontinuity. 
Certainly,  in  all  the  examples,  this  turning  point  becomes  sharper  with  increased 
compensator  order.  If  it  is  true  that  this  point  is  a  point  of  discontinuity,  it 
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immediately  follows  that  the  compensator  must  have  infinite  order.  The 
controller  (at  least  for  a  SISO  system)  is  nothing  more  than  a  ratio  of 
polynomials  in  the  Laplace  domain.  Therefore,  in  order  to  make  this 
discontinuous  turn,  an  infinite  number  of  polynomial  terms  are  required. 

Evidence  3:  Recall  that  two  of  the  necessary  conditions,  Equations  (3.10)  and 
(3.11),  are  Lyapunov  equations. 

AQ2  +  Q2At  +  BWB  l  =  0  (6.1) 

ATX  +  XA  +  C2TC2  =  0  (6.2) 


Consider  the  following  lemma: 

Lemma  6.1:  Suppose  X  >  0,  Z  >  0,  (>/Z,A)  is  detectable  and 
AtX  +  XA  +  Z  =  0 
Then  A  is  stable. 

Moreover,  X  >  0  iff  (V Z,A)  is  observable. 

Proof:  For  the  first  conclusion,  see  ([Won85, 283-284],  Theorem  12.2). 
For  the  second  conclusion,  see  ([Won85,58],  Theorem  3.1). 
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Consider  the  SISO  example.  It  was  discovered  that  in  every  case  examined,  both 
X  and  Q2  were  positive  definite.  Therefore,  by  Lemma  6.1  and  its  dual,  it 
follows  that  (A,BU)  is  controllable  and  (Cy,A)  is  observable.  In  the  SISO 
example,  (A,BU)  is  controllable  and  (Cy,A)  is  observable.  It  follows,  then,  that 
the  compensator  is  also  controllable  and  observable  (i.e.  it  has  a  minimal 
realization).  It  was  interesting  to  note,  however,  that  while  all  the  eigenvalues 
of  X  and  Q2  were  positive  and  nonzero,  usually  only  three  had  significant 
magnitude.  Most  of  the  rest  of  the  eigenvalues  were  very  close  to  zero.  This 
corroborates  with  the  earlier  observation  that  many  of  the  higher  order 
compensator  poles  had  zeros  almost  right  on  top  of  them  and  yet,  due  to  oo-norm 
bound,  the  pole/ zero  cancellations  could  not  actually  be  made.  If  it  is  true  that 
the  mixed  solution  is  indeed  a  minimal  realization  for  orders  up  to  infinity,  it 
would  follow  that  the  optimal  order  of  the  mixed  H^H,*  solution  is  infinity. 

6.2  Recommendations  for  Future  Research 

Obviously,  the  complete  proof  of  optimal  order  remains  to  be  shown.  This 
work  has  shed  some  light  on  this  subject,  but  it  will  require  further  research 
before  a  final  answer  can  be  given  to  the  question  of  optimal  order. 

This  research  dealt  completely  with  compensators  whose  order  is  greater  than 
the  order  of  the  plant.  As  mentioned  at  the  beginning  of  this  thesis,  the  more 
practical  research  would  be  in  the  area  of  reduced  order  compensators.  This  is 
a  much  more  difficult  area  of  study,  but  it  will  be  interesting  to  see  the  plot  of 
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a *  versus  compensator  order  completed  by  extending  it  below  the  order  of  the 
plant. 

One  issue  that  was  looked  at  (but  not  resolved)  is  the  problem  of  uniqueness 
of  the  solution.  Recall  that  one  of  the  necessary  conditions  in  the  general  mixed 
problem  is  a  Lyapunov  equation  with  no  constant  term.  The  only  way  to  get  a 
nonzero  solution  to  this  equation  is  by  requiring  its  "A"  matrix  to  be  neutrally 
stable.  However,  this  leads  to  an  infinite  number  of  solutions  of  varying  rank. 
It  is  not  clear  which  of  these  solutions  should  be  chosen  or  why  the  numerical 
solution  converges  to  one  solution  over  another.  In  terms  of  looking  at  the 
characteristics  of  higher  order  solutions,  nonuniqueness  of  the  solution  does  not 
cause  much  trouble.  All  the  solutions  that  were  obtained  are  valid  and  the  plots 
shown  in  Chapter  V  are  correct.  However,  when  talking  about  optimal  order, 
nonuniqueness  of  the  solution  is  a  problem.  If  there  is  a  family  of  solutions  to 
the  mixed  problem,  it  may  be  true  that  a  different  (full  order  or  smaller  increased 
order)  compensator  could  achieve  the  same  results  as  the  higher  order  controllers 
shown  here. 

Research  needs  to  continue  into  finding  a  closed-form  solution  to  the  mixed 
problem.  It  may  be  true  that  no  closed-form  solution  exists.  If  this  is  the  case, 
a  more  time  efficient  method  for  numerically  solving  the  problem  needs  to  be 
developed.  Some  of  the  higher  order  solutions  literally  took  one  to  two  weeks 
of  almost  continuous  run  time  for  a  single  solution.  Before  mixed  ^/H,*  can 
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become  a  useful  design  method  from  a  practical  standpoint,  further  research  into 
solution  algorithms  is  required. 

6.3  Summary 

In  this  thesis,  the  effects  of  using  higher  order  compensators  in  mixed 
have  been  investigated.  The  mixed  problem  in  general  has  been  discussed  and 
the  motivations  for  performing  mixed  optimization  have  been  presented. 

The  issue  of  why  higher  order  controllers  are  important  was  addressed  and  the 
optimal  orders  of  related  optimization  problems  were  shown.  Then,  the  general 
mixed  problem  was  developed  using  a  Lagrange  multiplier  technique.  The 
necessary  conditions  for  a  minimum  were  derived  and  examined.  Next,  due  to 
numerical  difficulties  in  solving  the  true  mixed  problem,  a  suboptimal  problem 
was  developed  and  new  necessary  conditions  given.  After  showing  some  key  full 
order  results,  theoretical  results  for  the  higher  order  case  were  presented.  In 
particular,  the  key  proofs  that  were  given  include:  the  global  minimum  2-norm 
is  unachievable  under  output  feedback  for  certain  levels  of  y  regardless  of 
compensator  order;  the  solution  to  the  mixed  problem  lies  on  the  boundary  of  the 
oo -norm  constraint  for  this  same  range  of  y’s;  and,  the  suboptimal  mixed 
problem  converges  to  the  optimal  in  the  limit  for  higher  order  compensators. 
Then,  numerical  SISO  and  MIMO  examples  were  examined.  It  was  seen  that 
higher  order  compensators  do  produce  a  lower  2-norm  and  they  are  better  able 
to  recover  to  the  H2  solution  than  the  full  order  controllers.  Finally,  based  on 
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the  results  from  the  numerical  examples,  it  was  conjectured  that  the  optimal  order 
of  the  mixed  solution  is  infinite. 
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Appendix:  FORTRAN  Source  Code  of 
DFP  Algorithm 


This  appendix  contains  the  FORTRAN  source  code  for  the  Davidon- 
Fletcher-Powell  numerical  algorithm  used  to  solve  the  mixed 
optimization  problem.  All  code  required  to  run  the  program  is  included  here 
except  for  the  Riccati  equation  solver  and  eigenvalue  solver  routines.  These 
routines  are  readily  available  as  public  domain  software.  The  calls  to  these 
routines,  which  will  need  to  be  modified  depending  on  the  software  used,  are 
in  the  subroutine  EVALUF  and  are  marked  by  bold  italics. 

Below  is  a  summary  of  the  routines  included: 


Name 

Description 

page 

DIMEN. INC 

Separate  file  containing  array  dimensions 

A-2 

DFP 

Main  program 

A-3 

INPDAT 

Inputs  system  data 

A-6 

INITGS 

Inputs  initial  guess  for  compensator 

A-10 

FINDAL 

Calculates  DFP  Alpha  star 

A-12 

UPDATH 

Updates  the  H  matrix  and  S  vector 

A-16 

CKSTOP 

Determines  if  solution  is  converged 

A-18 

WRITER 

Writes  ouput  data  to  RESET  file 

A-19 

EVALUF 

Evaluates  the  value  of  the  cost  function 

A-21 

EVDELF 

Evaluates  the  derivatives  of  Laplacian 

A-30 

INPUT.DAT 

Sample  input  file 

A-35 

C*  *  *  *  *  *  ★  *  *  ★  ★  *  *  *  *  *  *  ★  *  *  *  *  Hr  *  *  *  ★  *  *  *  ★  ★  ★  *  *  *  *  *  ★  ★  *  ★  ★  *  *  ★  *  *  ★  ★  it  it  *  ★  ★  ★  *  ★  *  ★  ★  *  *  *  *  * 

c 

C  THIS  FILE  SHOULD  BE  ' INCLUDE' D  IN  THE  DFP  PROGRAM.  IT  CONTAINS 

C  ALL  THE  ARRAY  SIZES  NEEDED  IN  THE  DFP  DECLARATIONS. 

C 

C  ISTATE  =  LENGTH  OF  SYSTEM  STATE  VECTOR 

C  ESTATE  =  LENGTH  OF  COMPENSATOR  STATE  VECTOR 

C 

C  NUMBD  =  LENGTH  OF  EXOGENOUS  INPUT  VECTOR  D 

C  NUMBW  =  LENGTH  OF  EXOGENOUS  INPUT  VECTOR  W 

C  NUMBU  =  LENGTH  OF  CONTROLLED  INPUT  VECTOR  U 
C 

C  NUMBE  =  LENGTH  OF  CONTROLLED  OUTPUT  VECTOR  E 

C  NUMBZ  =  LENGTH  OF  CONTROLLED  OUTPUT  VECTOR  Z 

C  NUMBY  =  LENGTH  OF  MEASURED  OUTPUT  VECTOR  Y 

C****************************************************************** 

c 

INTEGER  I STATE , ESTATE , NUMBD , NUMBW , NUMBU , NUMBE , NUMBZ , NUMBY , 
TILDIM, NUMB 
C 

PARAMETER  (ISTATE  =  3, 

ESTATE  =  18, 

.  NUMBD  =  1, 

NUMBW  =  1, 

NUMBU  =  1, 

NUMBE  =  1, 

.  NUMBZ  =  1, 

NUMBY  =  1 ) 

PARAMETER 

(TILDIM  =  ISTATE+ESTATE, 

NUMB  =  ( ESTATE *ESTATE)  +  ( NUMBU*ESTATE )  + 

( NUMBY*ESTATE ) ) 

C****************************************************************** 
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oooo  ooo  nan  non 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 
C  This  program  uses  the  Davidon-Fletcher-Powell  C 
C  numerical  optimization  algorithm  to  solve  the  C 
C  suboptimal  mixed  H2/Hinf/Entropy  optimization  C 
C  problem.  That  is,  determine  (AC,BC,CC)  that  C 
C  minimizes  the  cost  functional:  C 

c  c 

C  J  =  ( 1-AMU )*tr(Q2*CZ'*CZ)  +  AMU*tr ( QINF*CE ' *CE)  C 
C  C 
C  where  0.0<  AMU  <1.0  C 
C  C 
C  and  Q2  is  the  real,  symmetric,  positive  semidefinite  C 
C  solution  to  the  Lyapunov  equation:  C 
C  ATIL*Q2  +  Q2*ATIL'  +  BWTIL*BWTIL '  =0  C 
C  C 
C  and  such  that  the  Riccati  equation:  C 
C  ATIL*QINF  +  QINF*ATIL '  +  C 
C  GAM2 INV*QINF*CETIL ' *CETIL*QINF  +  BDTIL*BDTIL '  =0  C 
C  has  a  real,  symmetric,  positive  semidefinite  solution.  C 
C  C 


ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

IMPLICIT  REAL *8 ( A-H,0-Z ) 

c 

C  -  INCLUDE  the  file  containing  all  the  array  dimensions. 

C  NOTE:  The  file  DIMEN. INC  must  be  modified  every  time 
C  there  is  a  change  in  the  size  of  the  system  or  compensator, 

C  and  this  program  must  be  re-compiled. 

C 

INCLUDE  'DIMEN. INC' 

-  DFP  algorithm  variables 

COMMON /HMATRX/  H { NUMB , NUMB ) ,  S(NUMB),  DELF ( NUMB ) ,  DELOLF ( NUMB ) 
COMMON/FLAGS/  IFLAG2,  ICNT,  TOLCHK,  CHECKSTOP 
COMMON/ PARAM/  AMU,  OMAMU,  GAMMA,  GAM2INV 

-  Say  lello  and  set  the  counters 

PRINT  * , ' THE  GREAT  AND  MIGHTY  DFP' 

PRINT  * 

PRINT  *, 'Ahead  one-half  impulse  power,  Mr  Crusher.' 

I COUNT  =  0 
JCOUNT  =  0 
I STOP  =  0 

-  Open  the  input  and  output  files 

OPEN ( 1 , FILE= ' INPUT . DAT ' ) 

OPEN ( 2, FILE*' CHECK. DAT' ) 

OPEN ( 8, FILE* ' RESET.DAT ' ) 

OPEN ( 9, FILE=' OUTPUT. DAT' ) 

-  Input  the  system  matrices  and  program  parameters 

PRINT  *, 'CALLING  INPDAT ' 

CALL  INPDAT 
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-  Input  the  initial  guess  for  the  compensator 

PRINT  *, 'CALLING  INITGS' 

CALL  INITGS 

-  Evaluate  the  derivatives  of  the  Lagrangian 

PRINT  *,  'CALLING  EVDELF ' 

CALL  EVDELF 

-  Initialize  the  vector  S 

PRINT  *, 'CALCULATING  S' 

DO  10  1=1, NUMB 
SUM=0 . 0D0 
DO  20  J=1 , NUMB 

SUM=SUM-H ( I , J ) *DELF ( J ) 

20  CONTINUE 

S ( I ) =SUM 
10  CONTINUE 

-  Initialize  ALPHA  STAR 

ALSTAR  =  1.0D-08 

-  All  inputs  and  initializations  are  complete.  Begin  the 
iterations. 

PRINT  *, 'Engage. . . ' 

PRINT  * 

PRINT  * 

»>  This  is  the  return  point  for  the  iteration  «< 

30  CONTINUE 

-  Update  the  counters 

I COUNT= I COUNT+ 1 
J  COUNT* J COUNT+ 1 

-  Calculate  ALPHA  STAR,  the  step  size  in  the  S  direction  which 
minimizes  the  function 

PRINT  *, 'CALLING  FINDAL' 

CALL  FINDAL(ALSTAR, FSTAR) 

-  Save  current  derivatives  into  last-pass  derivatives 

DO  40  1*1, NUMB 

DELOLF ( I ) =DELF { I ) 

40  CONTINUE 

-  Evaluate  the  derivatives  of  the  Lagrangian 

CALL  EVDELF 
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-  Update  the  variables  H  and  S 

CALL  UPDATH ( ALSTAR ) 

-  Check  for  convergence  of  the  solution 

CALL  CKSTOP ( I STOP , FSTAR ) 

-  Write  updates  to  user  terminal  every  iteration 

WRITE ( * , 50 )  ESTATE f GAMMA , ALSTAR, FSTAR, 

+  I COUNT, AMU, CHECKSTOP 

WRITE (2,50) ESTATE , GAMMA, ALSTAR , FSTAR , 

+  I COUNT , AMU , CHECESTOP 

50  FORMAT ( / , I 2 , '  State  ( gam= ' , F7 . 4 , ' ) ' , E19 . 10, F17 . 8, 

+  '  —»>  '  ,  13  ,  /  ,  2X,  '  (mu=  ' ,  F7 . 5 ,  '  )  '  , 

+  '  CHECESTOP: ' ,E19. 8) 

-  Write  data  to  output  files  every  10  iterations 

IF ( JCOUNT . GE . 10 ) THEN 
WRITE ( 9 , * ) 

WRITE (9,*)  'OUT  OF  FINDAL ALSTAR , FSTAR, 

+  '  ITERATION' , ICOUNT 

CALL  EVALUF ( FTEMP ) 

WRITE (9,60)  (DELF( II) , 11=1, NUMB) 

60  FORMAT (4E20. 12) 

CALL  WRITER (ICOUNT) 

JCOUNT=0 
END  IF 

-  Check  for  completion  of  the  program  (quit  after  300  iterations 
if  no  solution  has  been  reached) 

IF ( ICOUNT . LT . 300 . AND . I STOP . NE . 1 )  GOTO  30 

-  Write  final  results  to  the  output  files 

WRITE (9,*) 

WRITE (9,*)  'FINAL  VALUES  OF  THE  DERIVATIVES: 

WRITE (9,60)  (DELF(II) ,11=1, NUMB) 

CALL  WRITER ( ICOUNT) 

PRINT*,  - - A' - aaaaaaaa*aaaa, - - - 

-  Clean  up  shop  and  go  home 

CLOSE ( 1 ) 

CLOSE ( 2 ) 

CLOSE ( 8 ) 

CLOSE ( 9 ) 

STOP 

END 
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CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCr’CCCCCCCCCC 


c 

c 

c 

c 

c 

c 

c 


INPDAT 

This  subroutine  reads  the  data  for  the  system  state- 
space  matrices  and  the  program  parameters  from  an 
input  file 


C 

c 

c 

c 

c 

c 

c 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  INPDAT 
IMPLICIT  REAL *8 { A-H, O-Z ) 

CHARACTER* 50  CHAR 


C 

C  -  INCLUDE  the  file  containing  all  the  array  dimensions. 

C  NOTE:  The  file  DIMEN. INC  must  be  modified  every  time 
C  there  is  a  change  in  the  size  of  the  system  or  compensator, 

C  and  this  subroutine  must  be  re-compiled. 

C 

INCLUDE  'DIMEN. INC' 


-  State  space  system  matrices 

COMMON/SYSTEM/ 

A( ISTATE, ISTATE) , 

BU ( ISTATE , NUMBU ) ,  BD ( ISTATE , NUMBD ) ,  BW( ISTATE, NUMBW) , 
CY(NUMBY, ISTATE) ,  CE ( NUMBE, ISTATE ) ,  CZ ( NUMBZ , ISTATE ) , 
DYD ( NUMBY , NUMBD ) ,  DYW ( NUMBY , NUMBW) ,  DEU ( NUMBE , NUMBU ) , 
DZU(NUMBZ, NUMBU) 


-  Optimization  parameters 

COMMON /PARAM/  AMU,  OMAMU,  GAMMA,  GAM2INV 

-  DFP  algorithm  variables 

COMMON /HMATRX/  H ( NUMB , NUMB ) ,  S(NUMB),  DELF ( NUMB ) ,  DELOLF ( NUMB ) 
COMMON /MATRIX/  XMATOL ( NUMB ) ,  TOLSER 
COMMON/FLAGS/  IFLAG2 ,  ICNT,  TOLCHK,  CHECKSTOP 


-  Format  statements 


500  FORMAT ( A50 ) 

510  FORMAT ( 817 ) 

520  FORMAT (2D11. 3) 

530  FORMAT ( 1D11 . 3 ) 

540  FORMAT (8E15. 5) 

541  FORMAT (8E15. 5) 

C 

C  -  All  inputs  that  are  read  from  the  input  file  are  echoed 
C  to  the  output  files. 

C 

C  -  Input/output  the  title 
C 

READ (1,500)  CHAR 
WRITE (9, 500)  CHAR 
WRITE (8, 500)  CHAR 
READ( 1,500)  CHAR 
WRITE (9, 500)  CHAR 
WRITE (8, 500)  CHAR 
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-  Tnput/output  the  dimensions  of  the  matrices 

NOTE:  Dimensions  are  not  really  input  here;  they  are  in  the 
INCLUDE  file.  They  appear  here  for  convenience. 

READ (1,500)  CHAR 
WRITE (8, 500)  CHAR 
READ (1,500)  CHAR 
>”,?TTE  (  8 , 500  )  CHAR 
READ (1,500)  CHAR 

WRITE (8,510)  ISTATE , ESTATE , NUMBU, NUMBY , NUMBD , NUMBE , NUMBW, NUMBZ 

-  Input /output  the  parameters  GAMMA  and  AMU 

READ (1,500)  CHAR 

WRITE (8, 500)  CHAR 

READ( 1,500)  CHAR 

WRITE (8, 500)  CHAR 

READ (1,520)  GAMMA, AMU 

WRITE (8, 520)  GAMMA, AMU 

GAM2INV  =  1 . 0D0/ ( GAMMA* GAMMA) 

OMAMU  =  1.0D0  -  AMU 

-  Input/output  the  tolerance  for  the  1-D  search  and  checkstop 

READ (1,500)  CHAR 

WRITE (8, 500)  CHAR 

READ (1,500)  CHAR 

WRITE (8, 500)  CHAR 

READ (1,520)  TOLSER , TOLCHK 

WRITE (8, 520)  TOLSER, TOLCHK 

-  Input/output  the  system  A  matrix 

READ( 1,500)  CHAR 
WRITE (8, 500)  CHAR 
READ (1,500)  CHAR 
WRITE (8, 500)  CHAR 
DO  10  1=1, ISTATE 

READ (1,540)  ( A ( I, J) ,J=1, ISTATE) 

WRITE (8,541)  (A(I,J) ,J=1, ISTATE) 

10  CONTINUE 

-  Input/output  the  system  BU  matrix 
Note:  Input  file  contains  BU  transpose 

READ( 1,500)  CHAR 
WRITE (8, 500)  CHAR 
READ( 1,500)  CHAR 
WRITE (8, 500)  CHAR 
DO  20  1=1, NUMBU 

READ (1,540)  (BU(J, I) ,J=1, ISTATE) 

WRITE (8, 541)  ( BU ( J, I ) ,J=1, ISTATE) 

20  CONTINUE 

-  Input/output  the  system  BD  matrix 
Note:  Input  file  contains  BD  transpose 

READ( 1,500)  CHAR 
WRITE (8, 500)  CHAR 
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READ (1,500)  CHAR 
WRITE (8, 500)  CHAR 
DO  30  1=1 , NUMBD 

READ (1,540)  (BD( J, I) , J=l, ISTATE) 

WRITE (8, 541)  (BD( J, I) ,J=1, ISTATE) 
30  CONTINUE 
C 

C  -  Input/output  the  system  BW  matrix 
C  Note:  Input  file  contains  BW  transpose 
c 

READ (1,500)  CHAR 
WRITE (8, 500)  CHAR 
READ( 1,500)  CHAR 
WRITE (8, 500)  CHAR 
DO  40  1=1 , NUMBW 

READ ( 1 , 540 )  (BW(J, I) ,J=1, ISTATE) 

WRITE (8, 541)  (BW( J, I) ,J=1, ISTATE) 
40  CONTINUE 
C 

C  -  Input/output  the  system  CY  matrix 
C 

READ( 1,500)  CHAR 
WRITE (8, 500)  CHAR 
READ( 1,500)  CHAR 
WRITE (8, 500)  CHAR 
DO  50  1=1 , NUMBY 

READ ( 1 , 540 )  (CY(I, J) ,J=1, ISTATE) 
WRITE (8, 541)  ( CY ( I , J ) ,J=1, ISTATE) 

50  CONTINUE 

C 

C  -  Input/output  the  system  CE  matrix 
C 

READ (1,5 00)  CHAR 
WRITE (8, 500)  CHAR 
READ (1,5 00)  CHAR 
WRITE (8, 500)  CHAR 
DO  60  1=1 , NUMBE 

READ (1,540)  (CE(I,J) ,J=1, ISTATE) 
WRITE (8,541)  ( CE ( I , J ) ,J=1, ISTATE) 

60  CONTINUE 
C 

C  -  Input/output  the  system  CZ  matrix 
C 

READ( 1,500)  CHAR 
WRITE (8, 500)  CHAR 
READ (1,500)  CHAR 
WRITE (8, 500)  CHAR 
DO  70  1=1 , NUMBZ 

READ ( 1 , 540 )  (CZ(I, J) ,J=1, ISTATE) 
WRITE (8, 541)  ( CZ( I, J) , J=l, ISTATE) 

70  CONTINUE 
C 

C  -  Input/output  the  system  DYD  matrix 
C 

READ (1,500)  CHAR 
WRITE (8, 500)  CHAR 
READ (1,500)  CHAR 
WRITE (8, 500)  CHAR 
DO  80  1=1, NUMBY 

READ { 1 , 540 )  (DYD(I, J) ,J=1, NUMBD) 
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WRITE(8, 541)  { DYD ( I , J ) , J=1,NUMBD) 

80  CONTINUE 
C 

C  -  Input/output  the  system  DYW  matrix 
C 

READ (1,500)  CHAR 
WRITE (8, 500)  CHAR 
READ( 1,500)  CHAR 
WRITE (8, 500)  CHAR 
DO  90  1=1 , NUMBY 

READ (1,540)  (DYW(I,J) ,J=1, NUMBW ) 
WRITE (8,541)  (DYW( I, J) ,J=1, NUMBW) 
90  CONTINUE 
C 

C  -  Input/output  the  system  DEU  matrix 
C 

READ (1,500)  CHAR 
WRITE (8, 500)  CHAR 
READ (1,500)  CHAR 
WRITE (8, 500)  CHAR 
DO  100  1=1 , NUMBE 

READ (1,540)  ( DEU ( I , J ) , J=1 , NUMBU ) 
WRITE (8, 541)  (DEU ( I , J ) ,J=1, NUMBU) 
100  CONTINUE 
C 

C  -  Input/output  the  system  DZU  matrix 
C 

READ( 1,500)  CHAR 
WRITE (8, 500)  CHAR 
READ (1,5 00)  CHAR 
WRITE (8, 500)  CHAR 
DO  110  1=1 , NUMBZ 

READ (1,540)  (DZU(I, J) , J=l, NUMBU) 
WRITE (8, 541)  (DZU( I, J) , J=l, NUMBU) 
110  CONTINUE 
C 

C  -  Initialize  the  H  matrix  to  the  identity 

c 

DO  120  1=1, NUMB 

DO  120  J=1 , NUMB 
H( I , J) =0 . 0D0 
IF ( I . EQ. J) H ( I , J ) =1 . 0D0 
120  CONTINUE 
RETURN 
END 


A-9 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


c  c 

C  INITGS  C 
c  c 
C  This  subroutine  reads  the  data  for  the  initial  guess  C 
C  of  the  state-space  compensator  matrices  from  an  C 
C  input  file  C 
C  C 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  INITGS 
IMPLICIT  REAL*8(A-H,0-Z) 

CHARACTER* 50  CHAR 
C 

C  -  INCLUDE  the  file  containing  all  the  array  dimensions. 

C  NOTE:  The  file  DIMEN. INC  must  be  modified  every  time 
C  there  is  a  change  in  the  size  of  the  system  or  compensator, 

C  and  this  subroutine  must  be  re-compiled, 

c 

INCLUDE  'DIMEN. INC' 

C 

C  -  Compensator  system  matrices 

c 

COMMON/COMP/  AC ( KSTATE , KSTATE ) ,  BC ( KSTATE , NUMBY ) , 

CC ( NUMBU , KSTATE ) 

C 

500  FORMAT (A50) 

510  FORMAT ( 4D19 . 11 ) 

C 

C  -  Input  the  AC  matrix 
C 

IDONE=0 
J2  =  0 

DO  5  WHILE (IDONE.EQ.0) 

J1  =  J2+1 
J2  =  Jl+3 

IF  (J2.GE. KSTATE)  THEN 
J2  =  KSTATE 
I DONE  »  1 
ENDIF 

READ( 1,500)  CHAR 
WRITE (9, 500)  CHAR 
READ (1,500)  CHAR 
WRITE (9, 500)  CHAR 
DO  10  1=1, KSTATE 

READ (1,5 10)  ( AC ( I , J ) , J=J1 , J2 ) 

WRITE (9,510)  ( AC ( I , J ) , J= J1  J2 ) 

10  CONTINUE 

5  CONTINUE 
C 

C  -  Input  the  BC  matrix  (NO  TRANSPOSE) 

C 

READ( 1,500)  CHAR 
WRITE (9, 500)  CHAR 
READ (1,5 00)  CHAR 
WRITE (9, 500)  CHAR 
DO  20  1=1, KSTATE 

READ (1,510)  (BC(I, J) ,J=1, NUMBY) 

WRITE (9,510)  (BC(I, J) ,J=1, NUMBY) 

20  CONTINUE 
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-  Input  the  CC  matrix 

IDONE  =  0 
J2  =  0 

DO  35  WHILE ( IDONE. EQ.O) 

J1  =  J2+1 
J2  =  Jl+3 

IF  (J2.GE. ESTATE)  THEN 
J2  =  ESTATE 
IDONE  =  1 
ENDIF 

READ (1,500)  CHAR 
WRITE (9, 500)  CHAR 
READ( 1,500)  CHAR 
WRITE (9, 500)  CHAR 
DO  30  1*1 , NUMBU 

READ (1,510)  (CC(I, J) , J=J1, J2) 
WRITE (9, 510)  (CC(I,J), J=J1, J2) 

30  CONTINUE 

35  CONTINUE 
RETURN 
END 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

FINDAL  C 

c 

The  subroutine  calculates  the  ALPHA  that  minimizes  C 

the  function  F(X  +  ALPHA*S)  C 

C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  FINDAL (ALSTAR, FSTAR) 

IMPLICIT  REAL*8 ( A-H, O-Z ) 

-  INCLUDE  the  file  containing  all  the  array  dimensions. 

NOTE:  The  file  DIMEN. INC  must  be  modified  every  time 

there  is  a  change  in  the  size  of  the  system  or  compensator, 
and  this  subroutine  must  be  re-compiled. 

INCLUDE  'DIMEN. INC' 

-  DFP  algorithm  variables 
COMMON/MATRIX/  XMATOL { NUMB ) ,  TOLSER 

COMMON /HMATRX/  H ( NUMB , NUMB ) ,  S ( NUMB ) ,  DELF (NUMB) ,  DELOLF ( NUMB ) 

-  Compensator  system  matrices 

COMMON/COMP/  AC (ESTATE, ESTATE ) ,  BC (ESTATE, NUMBY ) , 

+  CC(NUMBU, ESTATE) 


-  XMAT  is  a  vector  containing  the  matrices  AC,  BC,  and  CC 

DIMENSION  XMAT (NUMB) 

EQUIVALENCE  ( XMAT ( 1 ) , AC ( 1 , 1 ) ) 

EQUIVALENCE  ( XMAT (ESTATE*ESTATE+1 ) ,BC(1, 1) ) 

EQUIVALENCE  ( XMAT ( ESTATE*ESTATE+ESTATE*NUMBY+1 ) ,CC(1,1) ) 

-  Save  last-pass  XMAT  vector 

DO  10  1=1, NUMB 

XMATOL ( I ) =XMAT ( I ) 

10  CONTINUE 

-  Initialize  starting  ALPHA  to  last-pass  ALPHA  STAR 
SHIFT=ALSTAR 

-  Identify  an  initial  region 

II  =  0 
DO  30  1=1,3 
C  PRINT  *, 'I  =  ' , I 

ALPHA® ( 1-1 ) *SHIFT 
DO  40  J=1 , NUMB 

C  PRINT  *, 'XMATOL: ' ,J,XMATOL(J) 

XMAT ( J ) =XMATOL ( J ) + ALPHA*  S ( J ) 

C  PRINT  * , ' XMAT ( J )  =  ',J,XMAT(J) 

40  CONTINUE 

CALL  EVALUF(FF) 

C  PRINT  *, 'FUNCTION  =  '  ,  FF 
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c 

C  -  Identify  three  points  in  the  region 
C 

II  =  II  +  1 
IF  (II.EQ.l)  THEN 
FO=FF 
FI  =  FF 

ALPHA1  =  ALPHA 

C  PRINT  * , '  FI , ALPHA1 : ' , FI , ALPHA1 

ELSE  IF  (II.EQ.2)  THEN 
F2  *  FF 

ALPHA2  =  ALPHA 

C  PRINT  * , ' F2 , ALPHA2 : ' , F2  , ALPHA2 

ELSE 

F3  =  FF 

ALPHA3  =  ALPHA 

C  PRINT  * , ' F3 , ALPHA3 : ' , F3 , ALPHA3 

END  IF 
30  CONTINUE 
C 

C  -  Expand  upper  bound  until  the  minimizing  ALPHA  is  within  the 
C  region 
C 

DO  35  WHILE  (F3.LT.F2) 

F2  =  F3 

ALPHA2  =  ALPHA3 

ALPHA3  *  ALPHA3  *  2 . 0D00 

DO  36  1=1, NUMB 

XMAT(I)  =  XMATOL ( I )  +  ALPHAS  *  S(I) 

36  CONTINUE 

CALL  EVALUF ( F3 ) 

35  CONTINUE 
C 

C  -  Return  point  if  the  ALPHA  needs  refining 
C 

50  CONTINUE 

C  WRITE (2, 1000)  ALPHA1 , ALPHA2 , ALPHA3 , ALPHA4 , ALPHA5 , 

+  F1,F2,F3,F4,F5 

C1000  FORMAT (5E16.5,/,5F16.10,/) 

C 

C  -  Determine  which  region  contains  the  min 
C 

IF ( FI . LT . F2 ) THEN 

C  WRITE (2, 1001)  ALPHA1 , ALPHA2 , ALPHA3 , ALPHA4 , ALPHA5 , 

+  F1,F2,F3,F4,F5 

C1001  FORMAT ( ' 1  '  , 5E16 . 5 , / , 5F16 . 10 , / ) 

C 

C  -  The  min  definitely  lies  between  ALPHA1  and  ALPHA2,  shrink  the 
C  search  area 

C  ALPHA1 - HERE - ALPHA2  - ALPHAS 

C 

ALPHA3=ALPHA2 

ALPHA2  * ( ALPHA3-ALPHA1 ) / 2 . 0D0 
F3  -  F2 

DO  60  1=1, NUMB 

XMAT(I)  =  XMATOL ( I )  +  ALPHA2  *  S ( I ) 

60  CONTINUE 

CALL  EVALUF (F2) 

ELSE 


C  WRITE (2, 1002)  ALPHA1 , ALPHA2 , ALPHA3 , ALPHA4 , ALPHA5 , 

+  F1,F2,F3,F4,F5 

C1002  FORMAT ( ' 2  '  , 5E16 . 5 , / , 5F16 . 10 , / ) 

C 

C  -  The  min  lies  to  the  left  or  right  of  ALPHA2,  determine  which 
C  side 
C 

ALPHA4= { ALPHA2 -ALPHA1 )/2. 0D0+ALPHA1 
DO  90  1=1, NUMB 

XMAT ( I ) =XMATOL ( I ) +ALPHA4  *  S ( I ) 

90  CONTINUE 

CALL  EVALUF ( F4 ) 

IF ( F4 . LT . F2 ) THEN 
C 

C  -  The  min  lies  between  ALPHA1  and  ALPHA2 

C  ALPHA1 - HERE - ALPHA2  - ALPHA3 

C 

ALPHA3=ALPHA2 
ALPHA2=ALPHA4 
F3  =  F2 
F2  =  F4 

C  WRITE (2, 1003)  ALPHA1 , ALPHA2 , ALPHA3 , ALPHA4 , ALPHA5 , 

+  F1,F2,F3,F4,F5 

C1003  FORMAT ( ' 3  '  , 5E16 . 5 , / , 5F16 . 10 , / ) 

ELSE 

ALPHA5= ( ALPHA3-ALPHA2 ) /2 . QD0+ALPHA2 
DO  100  1=1, NUMB 

XMAT ( I ) =XMATOL ( I ) +ALPHA5  *S ( I ) 

100  CONTINUE 

CALL  EVALUF (F5) 

IF ( F5 . LT. F2 ) THEN 
C 

C  -  The  min  lies  between  ALPHA2  and  ALPHA3 

C  ALPHA1 - ALPHA2 - HERE - ALPHA3 

C 

ALPHA1=ALPHA2 
ALPHA2=ALPHA5 
FI  =  F2 
F2  =  F5 

C  WRITE (2, 1004)  ALPHA1 , ALPHA2 , ALPHA3 , ALPHA4 , ALPHA5 , 

+  F1,F2,F3,F4,F5 

C1004  FORMAT ( ' 4  ’ , 5E16 . 5 , / , 5F16 . 10, / ) 

ELSE 

C 

C  -  The  min  lies  between  ALPHA  4  and  ALPHA  5 

C  ALPHA1 - ALPHA4 - HERE - ALPHAS - ALPHA3 

C 

ALPHA1=ALPHA4 
ALPHA3= ALPHAS 
FI  =  F4 
F3  =  F5 

C  WRITE (2, 1005)  ALPHA1 , ALPHA2 , ALPHA3 , ALPHA4 , ALPHA5 , 

+  F1,F2,F3,F4,F5 

C1005  FORMAT ( ' 5  ' , 5E16 . 5 , / , 5F16 . 10 , / ) 

END  IF 
END  IF 
END  IF 
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-  Check  for  converence  on  ALPHA 

IF ( ABS ( (ALPHA2-ALPHA1) /ALPHA2) . GT . TOLSER ) GOTO  50 

-  Update  XMAT  using  new  ALPHA  STAR 

ALSTAR=ALPHA1 
I F ( F2  . LT .FI) THEN 

IF ( F3 . LT . F2 ) THEN 
ALSTAR=ALPHA3 
ELSE 

ALSTAR=ALPHA2 

ENDIF 

ELSE 

IF ( F3 . LT . FI ) THEN 
ALSTAR=ALPHA3 
ENDIF 
ENDIF 

DO  110  1=1, NUMB 

XMAT { I ) =XMATOL ( I ) +ALSTAR*  S ( I ) 

110  CONTINUE 

CALL  EVALUF ( FSTAR ) 

RETURN 

END 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

UPDATH  C 

c 

This  subroutine  updates  the  variables  H  and  S  C 

C 


SUBROUTINE  UPDATH (ALSTAR) 

IMPLICIT  REAL* 8 ( A-H, O-Z ) 

-  INCLUDE  the  file  containing  all  the  array  dimensions. 

NOTE:  The  file  DIMEN. INC  must  be  modified  every  time 

there  is  a  change  in  the  size  of  the  system  or  compensator, 
and  this  subroutine  must  be  re-complied. 

INCLUDE  'DIMEN. INC' 


-  DFP  algorithm  variables 

COMMON /HMATRX/  H ( NUMB , NUMB ) ,  S(NUMB),  DELF (NUMB) ,  DELOLF ( NUMB ) 
DIMENSION  Z(NUMB),  AM (NUMB, NUMB ) ,  AN ( NUMB , NUMB ) ,  HY (NUMB) 

-  Update  Z 

DO  10  1=1, NUMB 

Z ( I ) =DELF ( I ) -DELOLF  ( I ) 

10  CONTINUE 
TR2  Z=0 . 0D0 

DO  11  1=1, NUMB 

TR2Z=Z ( I ) *Z ( I ) 

11  CONTINUE 

I F ( TR2  Z . EQ . 0 . 0D0 ) THEN 
ISTOP=l 
RETURN 
END  IF 

-  Find  the  denominator  for  M 

SUM=0 . 0D0 
DO  20  1=1, NUMB 

SUM=SUM+S ( I ) *Z ( I ) 

20  CONTINUE 

FACM=ALSTAR/SUM 

-  Store  M 

DO  30  1=1, NUMB 

DO  30  J=1 , NUMB 

AM ( I , J) =S ( I ) *S ( J) *FACM 
30  CONTINUE 

-  Find  the  denominator  for  N 

SUM=0 . 0D0 
DO  40  1=1, NUMB 

DO  40  J=1 , NUMB 

SUM=SUM+Z (I)*H(I,J)*Z(J) 

40  CONTINUE 
FACN*-SUM 
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c 

C  -  Calculate  N 
C 

DO  50  1=1, NUMB 
SUM=0 . 0D0 
DO  60  J=1 , NUMB 

SUM=SUM+H (I,J)*Z(J) 

60  CONTINUE 

HY ( I ) =SUM 
50  CONTINUE 

DO  70  1=1, NUMB 

DO  70  J=1 , NUMB 

AN  ( I ,  J )  =HY  ( I )  *HY  ( J )  /FACN 
70  CONTINUE 
C 

C  -  Update  H 
C 

DO  80  1=1, NUMB 

DO  80  J=1 , NUMB 

H ( I ,  J ) =H ( I , J ) +AM ( I , J ) +AN ( I ,  J  ) 

80  CONTINUE 
C 

C  -  Update  S 
C 

DO  90  1=1, NUMB 
SUM=0 . 0D0 
DO  100  J=1 , NUMB 

SUM=SUM-H ( I , J ) *DELF ( J ) 

100  CONTINUE 

S( I)=SUM 
90  CONTINUE 
C 

C  -  Check  to  insure  H  direction  will  decrease  function 
C 

SUM=0 . 0D0 
DO  200  1=1, NUMB 

SUM=SUM+S ( I ) *DELF ( I ) 

200  CONTINUE 

I F ( SUM . GT . 0 . 0D0 ) THEN 
WRITE (*,*)  'HAD  TO  UPDATE  H' 

WRITE (9,*)  'HAD  TO  UPDATE  H' 

DO  210  1=1, NUMB 
DO  210  J=1 , NUMB 
H ( I , J ) =0 . 0D0 
IF ( I . EQ .  J ) H ( I , J ) =1 . 0D0 
210  CONTINUE 
C 

C  -  Update  S 
C 

DO  290  1=1, NUMB 
SUM=0 . 0D0 
DO  300  J=1 , NUMB 

SUM=SUM-H ( I , J ) *DELF{ J) 

300  CONTINUE 

S ( I ) =SUM 
290  CONTINUE 
END  IF 
RETURN 
END 
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c  c 

C  CKSTOP  C 

c  c 

C  This  subroutine  determines  if  the  solution  converged  C 
C  C 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  CKSTOP( ISTOP, FSTAR) 

IMPLICIT  REAL *8 { A-H, O-Z ) 

C 

C  -  INCLUDE  the  file  containing  all  the  array  dimensions. 

C  NOTE:  The  file  DIMEN. INC  must  be  modified  every  time 
C  there  is  a  change  in  the  size  of  the  system  or  compensator, 

C  and  this  subroutine  must  be  re-compiled, 

c 

INCLUDE  'DIMEN. INC' 

C 

C  -  DFP  algorithm  variables 
C 

COMMON /HMATRX/  H ( NUMB , NUMB ) ,  S(NUMB),  DELF ( NUMB ) ,  DELOLF ( NUMB ) 
COMMON/FLAGS/  I FLAG 2 ,  ICNT,  TOLCHK,  CHECKSTOP 
C 

C  -  Set  the  STOP  flag  to  zero 
C 

ISTOP=0 
TR2Z=0.0D0 
DO  21  1=1, NUMB 

TR2  Z= ( DELF ( I ) -DELOLF ( I ) ) * ( DELF ( I ) -DELOLF ( I ) ) 

21  CONTINUE 

IF ( TR2Z . EQ. 0 . 0D0 ) THEN 
ISTOP=l 
RETURN 
END  IF 
C 

C  -  Calculate  the  magnitude  of  the  residuals 
C 

SUM=0 . 0D0 
DO  10  1=1, NUMB 

DO  10  J=1 , NUMB 

SUM=SUM+DELF (I)*H(I,J)*DELF(J) 

10  CONTINUE 

CHECKSTOP  =  ABS ( SUM/FSTAR) 

IF (CHECKSTOP.LT. TOLCHK) ISTOP=l 

RETURN 

END 
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c 

WRITER  c 

c 

This  subroutine  writes  output  date  to  the  RESET  file  C 

C 


SUBROUTINE  WRITER 
IMPLICIT  REAL  *  8 ( A-H , O- Z ) 

-  INCLUDE  the  file  containing  all  the  array  dimensions. 

NOTE:  The  file  DIMEN. IND  must  be  modified  every  time 

there  is  a  change  in  the  size  of  the  system  or  compensator, 
and  this  subroutine  must  be  re-complied. 

INCLUDE  'DIMEN. INC' 


-  Compensator  system  matrices 

COMMON/COMP/  AC (ESTATE, ESTATE ) ,  BC ( ESTATE , NUMBY ) , 

+  CC(NUMBU, ESTATE) 

-  DFP  algorithm  variables 

COMMON /HMATRX/  H ( NUMB , NUMB ) ,  S(NUMB),  DELF ( NUMB ) ,  DELOLF ( NUMB ) 

-  Assorted  and  various  variables 


COMMON/OUTDAT/  TRACE,  TWONORM 


500 

FORMAT (A50) 

501 

FORMAT (  /  ,  '  THE 

AC 

MATRIX 

(COLUMNS  ',13,' 

502 

FORMAT (  / ,  '  THE 

BC 

MATRIX' ) 

503 

510 

FORMAT (  /  ,  '  THE 

FORMAT (4D19. 11) 

CC 

MATRIX 

(COLUMNS  ',13,' 

M3,')') 

M3,')') 


-  Output  the  AC  matrix 

IDONE  =  0 
J2  =  0 

DO  10  WHILE ( IDONE. EQ.O) 

J1  =  J2  +  1 
J2  =  Jl+3 

IF  (J2.GE. ESTATE)  THEN 
J2  =  ESTATE 
IDONE  *  1 
ENDIF 

WRITE (8,501)  J1,J2 
DO  20  1=1, ESTATE 

WRITE (8,510)  (AC ( I , J ) ,J=J1,J2) 
20  CONTINUE 

10  CONTINUE 


-  Output  the  BC  matrix  (NO  TRANSPOSE) 

WRITE (8, 502) 

DO  30  1=1, ESTATE 

WRITE (8,510)  (BC(I,J) , J=1 , NUMBY ) 
30  CONTINUE 
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c 

C  -  Output  the  CC  matrix 
C 

I DONE  =  0 
J2  =  0 

DO  40  WHILE (IDONE.EQ.O) 

J1  =  J2+1 
J2  =  Jl+3 

IF  (J2.GE. ESTATE)  THEN 
J2  =  ESTATE 
IDONE  =  1 
ENDIF 

WRITE (8,503)  J1 , J2 
DO  50  1  =  1 , NUMBU 

WRITE ( 8 , 510 )  (CC{I, J) , J=J1, J2) 
50  CONTINUE 

40  CONTINUE 
C 

WRITE (9,*)  'The  two  norm  is',TWONORM 
WRITE ( * »  * )  'The  two  norm  is',TWONORM 
RETURN 
END 
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cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

c 

EVALUF  C 

c 

This  subroutine  evaluates  the  value  of  the  cost  C 

function  C 

C 

ccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 


SUBROUTINE  EVALUF (FU) 
IMPLICIT  REAL *8 { A-H, 0-Z ) 


-  INCLUDE  the  file  containing  all  the  array  dimensions. 

NOTE:  The  file  DIMEN. INC  must  be  modified  every  time 

there  is  a  change  in  the  size  of  the  system  or  compensator, 
and  this  subroutine  must  be  re-compiled. 


INCLUDE  'DIMEN. INC' 


-  State  space  system  matrices 


COMMON/SYSTEM/ 

A( ISTATE, ISTATE) , 
BU ( I STATE , NUMBU ) , 
CY(NUMBY, ISTATE) , 
DYD ( NUMBY , NUMBD ) , 
D2U(NUMBZ, NUMBU) 


BD( ISTATE, NUMBD) , 
CE(NUMBE, ISTATE) , 
DYW ( NUMBY .NUMBW) , 


BW ( I STATE , NUMBW ) , 
CZ(NUMBZ, ISTATE) , 
DEU ( NUMBE , NUMBU ) , 


-  Tilde  matrices 


COMMON /QTWOTL/ 

QTWO 1 ( I STATE , I STATE ) , 
.  QTW021(KSTATE, ISTATE) 
COMMON /QINFTL/ 

QINF1 ( ISTATE , ISTATE ) , 
.  QINF21( KSTATE, ISTATE) 
COMMON /XTL/ 

XTL1 (ISTATE, ISTATE) , 
XTL2 1 ( KSTATE , ISTATE ) , 
COMMON /YTL/ 

YTL1( ISTATE, ISTATE) , 
YTL21 (KSTATE, ISTATE) , 


QTWO 12 ( ISTATE , KSTATE ) 
QTW02 ( KSTATE , KSTATE ) 

QINF12 ( ISTATE , KSTATE ) 
QINF2 ( KSTATE , KSTATE ) 

XTL12( ISTATE, KSTATE) , 
XTL2 ( KSTATE , KSTATE ) 

YTL12{ ISTATE, KSTATE) , 
YTL2 ( KSTATE , KSTATE ) 


-  Compensator  system  matrices 


COMMON/COMP/  AC (KSTATE, KSTATE ) ,  BC ( KSTATE , NUMBY ) , 
+  CC( NUMBU, KSTATE) 

-  Optimization  parameters 

COMMON 'PARAM/  AMU,  OMAMU,  GAMMA,  GAM2INV 


-  DFP  algorithm  variables 


COMMON /MATRIX/  XMATOL ( NUMB ) ,  TOLSER 
-  Riccati  solution  matrices 
COMMON/ RICINP/ 

F ( TILDIM, TILDIM ) ,  G ( TILD IM , TILDIM ) ,  H (TILDIM, TILDIM ) , 
X(TILDIM,TILDIM) 
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COMMON/RICSCR/ 

Z(4*TILDIM*TILDIM) ,  W(4*TILDIM*TILDIM) , 
ER(2*TILDIM) ,  El ( 2*TILDIM) , 

WORK ( 2 *TILDIM) ,  IND ( 2  *TILDIM) 

COMMON /RSOLN/ 

QTLINF(TILDIM, TILDIM) ,  QTLTWO ( TILDIM, TILDIM) , 
XTILDE (TILDIM, TILDIM) ,  YTILDE ( TILDIM, TILDIM) 


-  Assorted  and  various  variables 
COMMON/OUTDAT/  TRACE,  TWONORM 

COMMON/FLAGS/  IFLAG2,  ICNT,  TOLCHK,  CHECKSTOP 

COMMON /TILDES/  RTLTWO(TILDIM, TILDIM) ,  RTLINF ( TILDIM, TILDIM) , 
VTLTWO(TILDIM, TILDIM) ,  VTLINF ( TILDIM, TILDIM ) , 

ATIL(TILDIM, TILDIM)  , 

DTIL( TILDIM, TILDIM) , BWTIL ( TILD IM, TILDIM ) , 

CETIL(TILDIM, TILDIM) ,  CZTIL ( TILDIM, TILDIM) 

★  **★★★**•**********★***********★*★******★****★**★*★******★★***★*★ 

N  =  TILDIM 
NDIM  =  TILDIM 
ICNT  -  ICNT+1 
C 

C  -  Calculate  ATILDE 
C  ATIL=  (A  BU*CC] 

C  (BC*CY  AC] 

C 

DO  40  1=1 , I STATE 
DO  40  J=1 , ISTATE 

ATIL ( I , J )  =  A ( I , J ) 

40  CONTINUE 

DO  41  1=1, ESTATE 
DO  41  J=l, ESTATE 

ATIL ( 1+ ISTATE , J+ I STATE )  =  AC(I,J) 

41  CONTINUE 

DO  42  1=1, ISTATE 
DO  42  J=l, ESTATE 
SUM1  =  0.0D0 

DO  43  K=1 , NUMBU 

SUM1  =  SUM1+BU ( I , K ) *CC ( K,  J) 

43  CONTINUE 

ATIL ( I , J+ I STATE )  =  SUM1 

42  CONTINUE 

DO  44  1=1, ESTATE 
DO  44  J=l, ISTATE 
SUM1  =  0.0D0 

DO  45  K=1 , NUMBY 

SUM1  =  SUM1+BC ( I , K) *CY { K, J ) 

45  CONTINUE 

ATIL ( I+ISTATE , J )  =  SUM1 

44  CONTINUE 

C  WRITE (2,*) 

C  WRITE (2,*) 

C  WRITE] 2, *) 'ATIL' 

C  DO  46  1=1, TILDIM 

C  WRITE(2, 1000)  (ATIL(I, J) ,J=1, TILDIM) 

C  46  CONTINUE 
1000  FORMAT] 10E 15. 4) 
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-  Calculate  the  B  TILDES 

BDTIL  =  (  BD  J  BWTIL  =  [  BW  ] 

[BC*DYD]  ( BC+DYW] 

DO  70  1=1 , I STATE 
DO  70  J=1 , NUMBD 

BDTIL(I, J)=BD(I,  J) 

70  CONTINUE 

DO  72  1=1 , ISTATE 
DO  72  J=1 , NUMBW 

BWTIL ( I,J)=BW(I,J) 

72  CONTINUE 

DO  73  1=1, ESTATE 
DO  74  J=l, NUMBD 
SUM1  =  0.0D0 
DO  75  K=1 , NUMBY 

SUM1  =  SUM1+BC ( I ,  K) *DYD ( K, J ) 

75  CONTINUE 

BDTIL ( I+ISTATE , J )  =  SUM1 
74  CONTINUE 

DO  76  J=l, NUMBW 
SUM1  =  0.0D0 
DO  77  K=1 , NUMBY 

SUM1  =  SUM1+BC ( I , K) *DYW ( K, J ) 

77  CONTINUE 

BWTIL ( I+ISTATE, J)  =  SUM1 

76  CONTINUE 

73  CONTINUE 
WRITE ( 2 , * ) 

WRITE (2 , * ) 'BWTIL' 

DO  78  1=1 , TILDIM 

WRITE (2, 1000)  ( BWTIL ( I, J) ,J=1, NUMBW) 

78  CONTINUE 
WRITE (2, *) 

WRITE (2, *) 'BDTIL' 

DO  79  1=1, TILDIM 

WRITE (2, 1000)  ( BDTIL ( I, J) ,J=1, NUMBD) 

79  CONTINUE 

-  Calculate  the  C  TILDES 

CETIL  =  [CE  DEU*CC )  CZTIL  =  [CZ 

DO  80  1=1, ISTATE 
DO  81  J=1 , NUMBE 


CETIL{ J, I )  = 

CE( J, I ) 

81 

CONTINUE 

DO  82  J=1 , NUMBZ 

CZTIL( J, I )  = 

CZ(J,I) 

82 

CONTINUE 

80 

CONTINUE 

DO  83  1=1. ESTATE 
DO  84  J=l, NUMBE 
SUM1  =  0.0D0 
DO  85  K= 1 , NUMBU 

SUM1=  SUM1+DEU ( J, K) *CC ( K, I ) 
85  CONTINUE 

CETIL(J, I+ISTATE)  =  SUM1 
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84  CONTINUE 

DO  86  J=l, NUMBZ 
SUM1  =  0.0D0 
DO  87  K=1 / NUMBU 

SUM1  =  SUM1+DZU ( J,  K) *CC ( K,  I ) 

87  CONTINUE 
CZTIL( J, I+ISTATE)  =  SUM! 

86  CONTINUE 

83  CONTINUE 
WRITE (2, *) 

WRITE (2,*) ' CZTIL ' 

DO  88  1=1 , NUMBZ 

WRITE (2, 1000)  (CZTIL{I,J) ,  J=1 , TILDIM) 

88  CONTINUE 
WRITE (2 , * ) 

WRITE (2, *) 'CETIL' 

DO  89  1  =  1 , NUMBE 

WRITE (2, 1000)  ( CETIL ( I ,  J) ,  J=l, TILDIM) 

89  CONTINUE 

-  Calculate  the  R  TILDES 

RTLTWO=CZ  TILDE'  *  CZ  TILDE 
RTLINF=CE  TILDE'  *  CE  TILDE 

DO  160  1=1, TILDIM 
DO  160  J=l, TILDIM 
SUM1=0 • 0D0 
DO  170  K=l, NUMBZ 

SUM1=SUM1+CZTIL(K, I) *CZTIL(K, J) 

170  CONTINUE 

RTLTWO ( I , J ) =SUM1 
SUM1=0 . 0D0 
DO  180  K=l, NUMBE 

SUM1=SUM1+CETIL<K, I) *CETIL(K, J) 

180  CONTINUE 

RTLINF ( I , J) =SUM1 
160  CONTINUE 

C  WRITE (2,*) 

C  WRITE (2,*) 'RTLTWO' 

C  DO  181  1=1, TILDIM 

C  WRITE ( 2 , 1000 )  (RTLTWO(I,J) ,J=1, TILDIM) 

C  181  CONTINUE 
C  WRITE (2,*) 

C  WRITE (2, *) 'RTLINF' 

C  DO  182  1=1, TILDIM 

C  WRITE(2, 1000)  (RTLINFf I, J) ,J=1, TILDIM) 

C  182  CONTINUE 
C 

C  -  Calculate  the  V  TILDES 
C 

C  VTLTWO=BW  TILDE  *  BW  TILDE' 

C  VTLINF=BD  TILDE  *  BD  TILDE' 

C 

DO  190  1=1, TILDIM 
DO  190  J=l, TILDIM 
SUM1=0 . 0D0 
DO  200  K=1 , NUMBW 

SUM1=SUM1+BWTIL ( I , K) *BWTIL( J,K; 

200  CONTINUE 
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VTLTWO ( I , J ) =SUM1 
SUM1=0 . 0D0 
DO  210  K=1 , NUMBD 

SUM1=SUMI+BDTIL(I,K) *BDTIL( J,K) 

210  CONTINUE 

VTLINF ( I ,  J ) =SUM1 
190  CONTINUE 
WRITE (2,*) 

WRITE (2,*)' VTLTWO ' 

DO  211  1=1 , TILDIM 

WRITE(2, 1000)  (VTLTWO(I, J) ,  J=l, TILDIM) 

211  CONTINUE 
WRITE ( 2 , * ) 

WRITE ( 2 , * ) 'VTLINF' 

DO  212  1=1, TILDIM 

WRITE (2, 1000)  (VTLINF(I,J) ,J=1, TILDIM) 

212  CONTINUE 


-  Use  Riccati  solver 

(solving  F'X  +  XF  -  XGX  +  H  =  0  FOR  X) 

IFAIL=0 

Q2  SOLUTION 

F'X  +  XF  -  XGX  +  H  =0 

ATILDE*Q2  +  Q2*ATILDE '  -  X*0*X  +  BW*BW'  =0 

thus  F=ATILDE ' 

G=0 

H=BW*BW ' 

NOTE:  must  "start"  RICSOL  with  X=H 

DO  220  1=1, TILDIM 

DO  220  J=l, TILDIM 
F  ( I , J ) =ATIL ( J , I ) 

G ( I , J ) =0 . 0D0 
X ( I , J ) = VTLTWO ( I , J ) 

220  CONTINUE 

-  Call  Riccati  solver  for  solution  to  Q2 

PRINT  *, 'CALLING  RICSOL  1' 

CALL  RICCATI  SOLVER  (output:  X) 

PRINT  * , ' BACK  FROM  RICSOL  1' 

DO  230  1=1, TILDIM 

DO  230  J=l, TILDIM 

QTLTWO ( I , J ) =X ( I , J ) 

230  CONTINUE 
C  WRITE ( 2 , * ) 

C  WRITE (2,*) 'QTLTWO' 

C  DO  231  1=1, TILDIM 

C  WRITE(2, 1000)  (QTLTWO(I,J) ,J=1, TILDIM) 

C  231  CONTINUE 
C 

C  -  Check  for  stable  and  unique  solution 

c 

CALL  EIGENVALUE  SOLVER  (output:  ER-real  part  of  eig'B) 


A-25 


oooooooo  oo  oooo  ooooooooooo 


DO  301  1*1 , N 

C  WRITE ( 2 , * )  'Q2  EIGENV-REAL  (EIGINDEX) = ' , I , ' EIG* ' , ER< I ) 

IF(ER(I) .LE.-1E-20)  THEN 

C  PRINT*, 'Q2  PROBLEM  ( EIGINDEX) =', I EIG= ER( I ) 

IFAIL=1 
GO  TO  306 
ENDIF 
301  CONTINUE 


QINF  SOLUTION 

F'X  +  XF  -  XGX  +  H  =0 

ATILDE*QINF  +  QINF*ATILDE '  +  QINF*GAM2INV*CE ' *CE*QINF  +  BD*BD'=0 
thus  F=ATILDE ' 

G=-GAM2INV*CE ' *CE 
H=BD*BD ' 


NOTE:  F  is  the  same  as  for  Q2 
DO  240  1=1 , TILDIM 

DO  240  J=1,TILDIM 

G(I, J)=-GAM2INV  *  RTLINF ( I , J ) 

X(I, J)=VTLINF(I, J) 

240  CONTINUE 

-  Call  Riccati  solver  for  solution  to  QINF 

PRINT  *,  'CALLING  RICSOL  2' 

CALL  RICCATI  SOLVER  (output:  X) 

PRINT  * ,  '  BACK  FROM  RICSOL  2' 

DO  250  1*1, TILDIM 

DO  250  J=1 , TILDIM 

QTLINF (I,J)=X(I,J) 

250  CONTINUE 
WRITE ( 2 , * ) 

WRITE ( 2 , * ) 'QTLINF' 

DO  251  1=1, TILDIM 

WRITE(2, 1000)  (QTLINF(I, J) ,J*1, TILDIM) 

251  CONTINUE 

-  Check  for  stable  and  unique  solution 

CALL  EIGENVALUE  SOLVER  (output:  ER-real  part  of  eig'a) 

DO  302  1*1, N 

C  WRITE ( 2 , * )  'QINF  EIGENV-REAL  ( EIGINDEX )=', I ,' EIG= ', ER ( I ) 

IF ( ER( I ) . LE . -IE-20 )  THEN 

C  PRINT*, 'QINF  PROBLEM  ( EIGINDEX )=', I ,' EIG* ', ER( I ) 

IFAIL=1 
GO  TO  306 
ENDIF 
302  CONTINUE 


A-26 


no  o  o  nnnnnnnnn  oo  oooo  ooooooooooo 


X  LAGRANGE  MULTIPLIER  SOLUTION 


F'X  +  XF  -  XGX  +  H  =0 

ATILDE ' *X  +  X*ATILDE  -  X*0*X  +  { 1-AMU ) *CZ ' *CZ  =0 
thus  F=ATILDE 
G=0 

H=( 1-AMU) *CZ'*CZ 

I F ( I FLAG 2 . EQ . 1 )  THEN 
DO  260  1=1, TILDIM 

DO  260  J=1 , TILDIM 
F<I, J)=ATIL(I, J) 

G { I , J ) =0 . 0D0 
X  ( I ,  J ) =OMAMU*RTLTWO ( I , J) 

260  CONTINUE 

-  Call  Riccati  solver  for  solution  of  XTILDE 

PRINT  *,  'CALLING  RICSOL  3' 

CALL  RICCATI  SOLVER  (output:  X) 

PRINT  * , ' BACK  FROM  RICSOL  3' 

DO  270  1=1 , TILDIM 

DO  270  J=l, TILDIM 

XTILDE ( I , J ) =X ( I , J ) 

270  CONTINUE 

WRITE (2,*) 

WRITE (2, *) 'XTILDE' 

DO  271  1  =  1 ,  TILDIM 

WRITE (2, 1000)  <XTILDE(I, J) ,  J«l, TILDIM) 

271  CONTINUE 
END  IF 

-  Check  for  stable  and  unique  solution 

CALL  EIGENVALUE  SOLVER  (output:  ER-real  part  of  erg's) 

DO  303  1=1, N 

WRITE (2,*)  'XTILDE  EIGENV-REAL  ( EIG INDEX ) = ' , I , ' EIG= ' , ER ( I ) 
IF ( ER (I) .LE.-1E-10)  THEN 

PRINT*, 'XTILDE  PROBLEM  { EIGINDEX) EIG= ' , ER ( I ) 
IFAIL=1 
GO  TO  306 
ENDIF 
303  CONTINUE 


C  Y  LAGRANGE  MULTIPLIER  SOLUTION 
C 

C  F'X  +  XF  -  XGX  +  H  =0 

C  P ' *Y  +  Y*P  -  X*0*X  +  AMU*CE ' *CE  =0 
C  where  P=ATILDE+ ( GAM2 INV*QINF*CE ' *CE ) 

C  thus  F=P 

C  G=0 

C  H=AMU*CE ’ *CE 

C 
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DO  280  1=1 , TILDIM 

DO  280  J=1 , TILDIM 
SUM1=0 . 0D0 
G ( I , J) =0 . 0D0 
X { I , J ) =AMU*RTLINF ( I , J) 

DO  290  K=l, TILDIM 

SUM1=SUM1  +  QTLINF(I,K)*RTLINF(K, J) 

290  CONTINUE 

SUM1=SUM1*GAM2 INV  +  ATIL(I,J) 

F ( I , J) =SUM1 
280  CONTINUE 

-  Call  Riccati  solver  for  solution  to  YTILDE 

PRINT  *, 'CALLING  RICSOL  4' 

CALL  RICCATI  SOLVER  (output:  X) 

PRINT  * ,  '  BACK  FROM  RICSOL  4' 

DO  300  1  =  1 , TILDIM 

DO  300  J=l, TILDIM 

YTILDE ( I , J ) =X ( I , J ) 

300  CONTINUE 

WRITE (  2  ,  *  ) 

WRITE (2 , *) 'YTILDE' 

DO  299  1=1, TILDIM 

WRITE (2, 1000)  (YTILDE ( I, J) ,J=1, TILDIM) 

299  CONTINUE 

-  Check  for  stable  and  unique  solution 

CALL  EIGENVALUE  SOLVER  (output:  BR-real  part  of  eig’B) 

DO  304  1=1, N 

WRITE ( 2 , * )  'YTILDE  EIGENV-REAL  ( EIG INDEX ) = ' , I , ' EIG= ' , ER( I ) 
IF(ER(I) .LE.-1E-10)  THEN 

PRINT*,  'YTILDE  PROBLEM  ( EIGINDEX)  =  ' , I ,  'EIG= ' , ER ( I ) 
IFAIL=1 
GO  TO  306 
END  IF 

304  CONTINUE 

CALL  EIGENVALUE  SOLVER  (output:  ER-real  part  of  eig'a) 

DO  305  1=1, N 

WRITE ( 2 , * )  'Ast  EIGENV-REAL  ( EIGINDEX )=', I ,' EIG= ', ER ( I ) 
IF(ER(I) .GE.-1E-8)  THEN 

PRINT*, 'Ast  PROBLEM  ( EIGINDEX )=', I ,' EIG= ', ER ( I ) 

IFAIL=1 
GO  TO  306 
END  IF 

305  CONTINUE 


-  If  acceptable  solutions  were  found  to  the  Riccati  and  Lyapunov 
equations,  use  these  solutions  to  calculate  the  function  value. 
Otherwise,  no  acceptable  solution  exists.  Set  the  function 
value  to  "very  big.” 


306 


C 


IF ( IFAIL . EQ. 1 ) THEN 
F  u  -  i  E  t  1 6 
ELSE 
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-  Calculate  the  cost  function 

PRINT  *, 'CALCULATING  THE  COST  FUNCTION' 

TRACE=0 . 0D0 

SUM=0 . 0D0 

DO  390  1  =  1 , TILDIM 

DO  390  J=1,TILDIM 

SUM=SUM+ ( 1 . ODO-AMU ) *QTLTWO ( I , J) *RTLTWO ( J , I ) 
+AMU*QTLINF ( I , J ) *RTLINF ( J ,  I ) 

390  CONTINUE 

FU=TRACE+SUM 
WRITE (2,*) 

WRITE (2, 1001)  FU 
1001  FORMAT (F10. 5) 

-  Calculate  the  2-norm  of  Tzw 

SUM=0 . 0D0 

DO  400  1=1, TILDIM 

DO  400  J=l, TILDIM 

SUM=SUM+QTLTWO ( I , J ) *RTLTWO ( J,  I ) 

400  CONTINUE 

TWONORM  =  DSQRT ( SUM) 

END  IF 
RETURN 
END 
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c  c 

C  EVDELF  C 

c  c 

C  This  subroutine  evaluates  derivatives  of  the  Laplacian  C 

C  with  respect  to  the  variable  matrices  AC,  BC,  and  CC  C 

C  C 


cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc 

SUBROUTINE  EVDELF 
IMPLICIT  REAL *8 (A-H,0-Z ) 

C 

C  -  INCLUDE  the  file  containing  all  the  array  dimensions. 

C  NOTE:  The  file  DIMEN. INC  must  be  modified  every  time 
C  there  is  a  change  in  the  size  of  the  system  or  compensator, 

C  and  this  subroutine  must  be  re-compiled. 

C 


c 

c 

c 


c 

c 

c 


c 

c 

c 


INCLUDE  'DIMEN. INC' 

DIMENSION  ACDER (ESTATE, ESTATE ) ,  BCDER ( ESTATE , NUMBY ) , 
CCDER ( NUMBU , ESTATE ) 

State  space  system  matrices 


COMMON /SYSTEM/ 

A ( ISTATE , ISTATE ) , 
BU( I STATE, NUMBU) , 
CY( NUMBY, ISTATE) , 
DYD ( NUMBY, NUMBD ) , 
DZU(NUMBZ, NUMBU) 


BD( ISTATE, NUMBD) , 
CE(NUMBE, ISTATE) , 
DYW( NUMBY, NUMBW) , 


BW( ISTATE, NUMBW) , 
CZ(NUMBZ, ISTATE) , 
DEU(NUMBE, NUMBU) , 


Tilde  matrices 


COMMON /QTWOTL/ 

QTWOl (ISTATE, ISTATE) ,  QTW012 ( ISTATE , ESTATE ) , 
QTW021 (ESTATE, ISTATE) ,  QTW02 ( ESTATE , ESTATE ) 
COMMON /QINFTL/ 

QINF1 (ISTATE, ISTATE) ,  QINF12 ( ISTATE , ESTATE ) , 
QINF21 (ESTATE, ISTATE) ,  Q I NF2 ( ESTATE , ESTATE ) 
COMMON /XTL/ 

XTL1 (ISTATE, ISTATE) ,  XTL12 ( ISTATE , ESTATE ) , 
XTL21 (ESTATE, ISTATE) ,  XTL2 (ESTATE, ESTATE) 
COMMON /YTL/ 

YTL1 (ISTATE, ISTATE) ,  YTL12 ( ISTATE , ESTATE ) , 
YTL21 (ESTATE, ISTATE) ,  YTL2 ( ESTATE , ESTATE ) 


Compensator  system  matrices 


COMMON/COMP/  AC (ESTATE, ESTATE ) ,  BC ( ESTATE, NUMBY ) , 
+  CC (NUMBU, ESTATE) 

C 

C  -  Optimization  parameters 
C 


COMMON/PARAM/  AMU,  OMAMU,  GAMMA,  GAM2INV 
C 

C  -  DFP  algorithm  variables 
C 


COMMON /HMATRX/  H ( NUMB , NUMB ) ,  S(NUMB),  DELF ( NUMB ) ,  DELOLF ( NUMB ) 
EQUIVALENCE  ( DELF ( 1 ) , ACDER ( 1 , 1 ) ) 

EQUIVALENCE  ( DELF ( ESTATE *ESTATE+U , BCDER (1,1) ) 

EQUIVALENCE  ( DELF ( ESTATE* ESTATE +ESTATE* NUMB Y+l ) ,CCDER(1,1) ) 
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-  Riccati  solution  matrices 

COMMON/RSOLN/ 

QTLINF ( TILDIM, TILDIM) ,  QTLTWO ( TILDIM, TILDIM) , 
.  XTILDE ( TILDIM, TILDIM) ,  YTILDE (TILDIM, TILDIM) 

-  Assorted  and  various  variables 

COMMON/FLAGS/  IFLAG2,  ICNT,  TOLCHK,  CHECKSTOP 


I FLAG 2  =  1 

PRINT  *, 'CALLING  EVALUF  FROM  EVDELF' 

CALL  EVALUF (FU) 

PRINT  * , ' BACK  FROM  EVALUF' 

I FLAG  =  0 

-  Partition  the  tilde  matrices  (XTILDE,  YTILDE,  QTLTWO,  QTLINF) 

Partition  the  1,1  components 

DO  10  1=1 , I STATE 
DO  10  J=1 , ISTATE 

XTL1 ( I ,  J )  =  XTILDE ( I, J) 

YTL1 ( I , J )  =  YTILDE ( I , J ) 

QTWOl ( I , J )  =  QTLTWO ( I, J) 

QINF1 ( I , J )  =  QTLINF ( I , J ) 

10  CONTINUE 

Partition  the  1,2  components 

DO  11  1=1, ISTATE 
DO  11  J=1 , KSTATE 

XTL12 ( I , J )  =  XTILDE (I, J+ISTATE) 

YTL12 ( I , J )  =  YTILDE (I, J+ISTATE) 

QTWOl 2 ( I , J )  =  QTLTWO ( I, J+ISTATE) 

QINF12 ( I , J )  =  QTLINF( I, J+ISTATE) 

11  CONTINUE 

Partition  the  components 

DO  12  1=1, KST 
DO  12  J=1 , 1ST 

XTL2 1 ( I , J )  XTILDE ( 1+ ISTATE , J ) 

YTL2 1 ( I , J )  =  YTILDE ( I+ISTATE , J) 

QTW02 1 ( I , J )  =  QTLTWO ( I+ISTATE, J) 

QINF21 ( I , J)  =  QTLINF ( I+ISTATE, J) 

12  CONTINUE 

Partition  the  2,2  components 

DO  13  1=1, KSTATE 
DO  13  J=l, KSTATE 

XTL2 ( I , J )  =  XTILDE( I+ISTATE, J+ISTATE) 

YTL2 ( I , J )  =  YTILDE< I+ISTATE, J+ISTATE) 

QTW02 ( I , J)  =  QTLTWO ( I+ISTATE, J+ISTATE) 

QINF2 ( I , J)  =  QTLINF ( I+ISTATE, J+ISTATE) 

13  CONTINUE 
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-  Find  the  derivative  wrt  AC 

ICOUNT=0 

DO  20  1=1 , KSTATE 

DO  20  J=l, KSTATE 
ICOUNT=ICOUNT+l 
SUM=0 . 0D0 
DO  25  K=1 , ISTATE 

SUM=SUM+XTL12 (K, I) *QTW021( J,K) 

SUM=SUM+XTL2 1 ( I , K ) *QTW012 ( K, J ) 
SUM=SUM+YTL12 ( K,  I ) *QINF21 ( J,K) 

SUM=SUM+YTL2 1  ( I ,  K ) *QINF12 ( K, J ) 

25  CONTINUE 

DO  30  K=l, KSTATE 

SUM=SUM+XTL2 (K, I) *QTW02 ( J,K) 

SUM=SUM+XTL2 ( I , K) *QTW02 ( K, J ) 

SUM=SUM+YTL2 ( I , K ) *QINF2 (K, J) 

SUM=SUM+YTL2 (K, I ) *QINF2 ( J,K) 

30  CONTINUE 

ACDER ( I , J) =SUM 
20  CONTINUE 

-  Find  the  derivative  wrt  BC 

DO  40  1=1, KSTATE 
DO  40  J=1 , NUMBY 

ICOUNT=ICOUNT+l 
SUM=0.0D0 
DO  45  K=l, ISTATE 
DO  45  L=l, ISTATE 

SUM=SUM+XTL 1 2 ( K , I ) *QTW01 (L,K)*CY(J,L) 
SUM=SUM+XTL2 1 ( I , K ) *QTW01 ( K, L ) *CY ( J , L) 
SUM=SUM+YTL12 (K, I) *QINF1(L,K) *CY( J,L) 
SUM=SUM+YTL2 1 ( I , K ) *QINF1 ( K, L ) *CY ( J, L ) 

45  CONTINUE 

DO  50  K=l, KSTATE 
DO  50  L=l, ISTATE 

SUM=SUM+XTL2 ( I , K) *QTW02 1(K,L)*CY(J,L) 
SUM=SUM+XTL2 ( K, I ) *QTW012 (L,K) *CY ( J , L ) 
SUM=SUM+YTI,2  ( I ,  K  )  *QINF2 1  (K,L)  *CY  ( J ,  L ) 
SUM=SUM+YTL2 ( K, I ) *QINF12 (L,K) *CY( J,L) 

50  CONTINUE 

DO  60  K=l, ISTATE 
DO  60  L=1 , NUMBW 

SUM=SUM+XTL12 (K, I) *BW(K,L)*DYW( J,L) 
SUM=SUM+XTL21( I,K) *DYW( J,L) *BW(K,L) 

60  CONTINUE 

DO  70  K=l, ISTATE 
DO  70  L=1,NUMBD 

SUM=SUM+YTL12 (K, I ) *BD ( K , L ) *DYD ( J, L ) 
SUM=SUM+YTL2 1 { I , K) *DYD ( J , L ) *BD ( K, L ) 

70  CONTINUE 

DO  80  K=l, KSTATE 
DO  80  L=l, NUMBY 
DO  80  M=l, NUMBW 

SUM=SUM+XTL2 ( I , K ) *BC ( K, L ) *DYW(L,M) *DYW( J,M) 
SUM=SUM+XTL2 (K, I ) *BC { K, L ) *DYW ( L , M ) *DYW ( J, M) 
80  CONTINUE 
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DO  90  K=1 , KSTATE 
DO  90  L=1 , NUMBY 
DO  90  M=1 , NUMBD 

SUM=SUM+YTL2 (I,K)*BC(K,L) *DYD ( L, M) *DYD( J,M) 
SUM=SUM+YTL2 (K, I ) *BC ( K, L) *DYD (L,  M) *DYD ( J,M) 

90  CONTINUE 

BCDER(I, J)=SUM 
40  CONTINUE 

-  Find  the  derivative  wrt  CC 

DO  100  1=1 , NUMBU 
DO  100  J=l, KSTATE 
ICOUNT=ICOUNT+l 
SUM=0 . 0D0 
DO  110  K=1 , I STATE 
DO  110  L=1 , I STATE 

SUM=SUM+BU ( K ,  I )  *  XTL1 ( L , K ) *QTW02 1 ( J , L ) 

SUM=SUM+BU(K, I) *XTL1(K,L)*QTW012(L, J) 

SUM=SUM+BU(K, I ) *YTL1 ( L , K ) *QINF2 1 ( J, L ) 

SUM=SUM+BU ( K, I ) *YTL1 { K,  L) *QINF12 ( L, J ) 

110  CONTINUE 

DO  111  K=1 , I STATE 
DO  111  L=l, KSTATE 

SUM=SUM+BU (K, I ) *XTL12 (K,L) *QTW02 (L, J) 

SUM=SUM+BU(K, I) *XTL21(L,K) *QTW02 ( J,L) 

SUM=SUM+BU(K, I ) *YTL12 (K,L) *QINF2(L, J) 

SUM=SUM+BU { K, I ) * YTL2 1 ( L , K ) *QINF2 ( J , L ) 

111  CONTINUE 

DO  120  K=1 , NUMBE 
DO  120  L=1 , ISTATE 

DO  121  K=l, ISTATE 
DO  121  N=l, ISTATE 
SUM=SUM+GAM2 INV* 

DEU ( K, I ) *CE ( K, L ) *QINF1 ( L , M ) *YTL1 (M, N) *QINF12 (N, J ) 
SUM=SUM+GAM2 INV* 

DEU (K, I)*CE(K,L) *QINF1 (M, L ) *YTL1 ( N, M) *QINF21 ( J, N ) 

121  CONTINUE 

DO  122  M=l, ISTATE 
DO  122  N=l, KSTATE 
SUM=SUM+GAM2 INV* 

DEU (K,I)*CE(K,L)*QINF1(M,L) *YTL2 1 ( N, M) *QINF2 { J, N ) 
SUM=SUM+GAM2 INV* 

DEU(K, I) *CE(K,L) *QINF1(L,M) *YTL12(M,N) *QINF2 (N, J) 

122  CONTINUE 

DO  123  M=l, KSTATE 
DO  123  N=l, ISTATE 
SUM*SUM+GAM2 INV* 

DEU { K, I ) *CE ( K, L ) *QINF2 1 ( M , L ) *YTL12 ( N, M) *QINF2 1 ( J, N ) 
SUM=SUM+GAM2 INV* 

DEU (K,I)*CE(K,L)*QINF12(L,M) *YTL2 1 ( M, N ) *QINF12 ( N, J) 

123  CONTINUE 

DO  124  M=l, KSTATE 
DO  124  N=l, KSTATE 
SUM=SUM+GAM2INV* 

DEU (K, I)*CE(K,L)*QINF12(L,M) *YTL2 (M, N ) *QINF2 ( N , J ) 
SUM=SUM+GAM2 INV* 

DEU{K, I) *CE(K,L) *QINF21 (M,L)*YTL2 (N,M)*QINF2( J,N) 

124  CONTINUE 
120  CONTINUE 
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DO  130  K=1 , NUMBE 
DO  130  L=1 , NUMBU 
DO  130  M=l, KSTATE 
DO  131  N=1 , ISTATE 
DO  131  0=1, ISTATE 

SUM=SUM+GAM2 INV*DEU ( K, I ) *DEU ( K, L ) * 

CC(L,M) *QINF21(M,N) * YTL1 ( N,  O) *QINF12 (0,J) 
SUM=SUM+GAM2INV*DEU{K, I ) *DEU(K,L) * 

CC ( L , M ) *QINF12 (N,M) * YTL1 ( O , N ) *QINF2 1 ( J , O ) 

131  CONTINUE 

DO  132  N=l, ISTATE 
DO  132  0=1, KSTATE 

SUM=SUM+GAM2INV*DEU(K, I) *DEU(K,L) * 

CC ( L, M) *QINF2 1 ( M, N ) * YTL12 ( N, O ) *QINF2 (O, J ) 
SUM=SUM+GAM2 INV*DEU { K, I ) *DEU ( K, L ) * 

.  CC(L,M) *QINF12(N,M) * YTL2 1 ( O , N ) *QINF2 ( J , O ) 

132  CONTINUE 

DO  133  N=l, KSTATE 
DO  133  0=1, ISTATE 

SUM=SUM+GAM2 INV*DEU (K, I ) *DEU{K,L) * 
CC(L,M)*QINF2(M,N) * YTL2 1 ( N , O ) *QINF12 (O, J) 
SUM=SUM+GAM2INV*DEU(K, I) *DEU(K,L) * 

CC ( L , M ) *QINF2 ( N , M ) * YTL12 ( O , N ) *QINF2 1 ( J , O ) 

133  CONTINUE 

DO  134  N=l, KSTATE 
DO  134  0=1, KSTATE 

SUM=SUM+GAM2INV*DEU(K, I) *DEU(K,L) * 

CC ( L , M ) *QINF2  <M,N) *YTL2<N,0) *QINF2 (O, J) 
SUM=SUM+GAM2 INV*DEU (K, I ) *DEU ( K, L ) * 

CC(L,M) *QINF2 (N,M) *YTL2 (0,N) *QINF2 ( J,0) 

134  CONTINUE 

130  CONTINUE 

DO  140  K=1 , NUMBZ 
DO  140  L=l, ISTATE 

SUM=SUM+ ( 1 . ODO-AMU )*DZU(K,I)*CZ(K,L) *QTW012 ( L , J ) 
SUM=SUM+( 1.0D0-AMU) *DZU(K, I) *CZ(K,L) *QTW021( J,L) 
140  CONTINUE 

DO  150  K=l, NUMBZ 
DO  150  L=l, NUMBU 
DO  150  M=l, KSTATE 

SUM=SUM+ ( 1 . ODO-AMU ) 

*DZU(K,I)*DZU(K,L)*CC(L,M) *QTW02 ( J,M) 
SUM=SUM+ ( 1 . ODO-AMU ) 

*DZU(K, I) *DZU(K,L) *CC(L,M) *QTW02 (M, J) 

150  CONTINUE 

DO  160  K=l, NUMBE 

DO  160  L=l, ISTATE 

SUM=SUM+AMU*DEU(K, I)*CE(K,L)*QINF12(L,J) 
SUM=SUM+AMU*DEU ( K, I ) *CE ( K , L ) *QINF2 1 ( J, L ) 

160  CONTINUE 

DO  170  K=l, NUMBE 
DO  170  L=l, NUMBU 
DO  170  M=l, KSTATE 

SUM=SUM+AMU*DEU(K, I )*DEU(K,L) *CC(L,M) *QINF2 ( J,M) 
SUM=SUM+AMU*DEU ( K, I ) *DEU ( K , L ) *CC ( L , M ) *QINF2 (M, J) 
170  CONTINUE 

CCDER ( I , J ) =SUM 
100  CONTINUE 
RETURN 
END 
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THIS  IS  THE  INPUT  FILE  FOR  THE  DIRECT  METHOD 
THE  SISO  MIX  8  STATE  COMPENSATOR  GAMMA  2.5 


THE  DIMENSIONS  I STATE , ESTATE , NU , N Y , ND , NE , NW , N Z 

38111111 

THE  PARAMETERS  GAMMA  AND  MU  (2D11.6) 

0.250D+01  0 . 100D-00 

THE  TOLERANCES  OF:  1-D  SEARCH,  CHECKSTOP  (2D11.6) 

0 . 100D-03  0 . 100D-07 

THE  A  MATRIX  (8F8.4) 

-0 . 39080E+00  -0 . 45650E+00  0.12657E+01 

0 . 14453E+01  -0 . 10491E+01  -0.12077E+01 

-0 . 12880E+00  0 . 67440E+00  0.10324E+01 

THE  BU  MATRIX  AS  BU  TRANSPOSE 

-0 . 427  50E+00  -0.44700E+00  -0.91720E+00 

THE  BD  MATRIX  AS  BD  TRANSPOSE 

0 . 48800E-01  0 . 36080E+00  0.35640E+00 

THE  BW  MATRIX  AS  BW  TRANSPOSE 

0 . 14077E+01  0 . 97230E+00  -0.16050E+01 

THE  CY  MATRIX 

-0 . 15567E+01  -0 . 19432E+01  -0.91400E-01 

THE  CE  MATRIX 

0 . 94200E+00  0 . 14400E-01  0.11870E+00 

THE  CZ  MATRIX 

-0 . 45000E-01  0 . 36060E+00  0.18972E+01 

THE  DYD  MATRIX 
0. 51850E+00 

THE  DYW  MATRIX 
0 . 38990E+00 

THE  DEU  MATRIX 
0. 13575E+01 

THE  DZU  MATRIX 
0. 57810E+00 


THE  AC  MATRIX  (COLUMNS  1  -  4) 

-0.42464048036D+01  -0. 4 193865 5996D+01 
0. 19735258928D+00  -0 . 13358453000D+01 
-0. 67612139041D+01  -0 . 52273410600D+01 
0.93570728936D+01  0. 84381238744D+01 

0. 34 126422 373D+01  0 . 30774899656D+01 

0.40025571618D+01  0 . 36094699 197D+0 1 

0. 56863512443D+00  0 . 5 1279002 246D+00 

-0. 1757029 12 39D+01  -0 . 15844730042D+01 


-0. 20982704551D+01  0 . j8597033962D+01 

-0.4882  5701 128D+01  0 . 4428003193  "'0+01 

-0. 62417707662D+01  0 . 84680779691D+01 

0.76522396137D+01  -0 . 22312 1159 10D+02 
0. 27908680859D+01  -0 . 69399026978D+01 
0. 32733021125D+01  -0 . 92716881072D+00 
0.46503134842D+00  -0. 20016352363D+02 
-0. 14369031873D+01  0 . 91291800491D+01 


THE  AC  MATRIX  (COLUMNS  5  -  8) 

0. 36629700415D+01  -0 . 7 62 2 2 70701 1D+00  -0 . 27 165890667D+01  -0 . 19927768328D+00 
0. 42023029687D+01  -0 . 87445680506D+00  -0 . 3 1 1 65775779D+01  -0 . 22861917803D+00 
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0.80364506611D+01 
-0 . 13945 6353 79D+02 
-0.67197475912D+01 
-0.35055010102D+01 
-0. 10543648466D+02 
0. 32035395034D+01 

THE  BC  MATRIX 
-0. 16294768664D+01 
0. 18976447418D+00 
-0. 2392362063 6D +01 
0.36761772931D+01 
0. 13407481212D+01 
0. 15725120366D+01 
0.22340357463D+00 
-0. 69029606181 D+00 


-0. 1672304215 4D+01 
0. 5134115515 6D+01 
0 . 9563483536 1D+00 
0. 11063088353D+01 
0.20770424268D+01 
-0.345 1245045 4D+01 


-0 . 596011808820+01 
0 . 84823843944D+01 
0 . 9 118822 002 5D+01 
0 . 93504776975D+01 
-0. 11037259724D+02 
0. 31119395199D+01 


0 . 8 548709883 3D+01 


0 . 72 133650879D+01 


-0.43720949159D+00 
-0.26400498205D+01 
0. 13756622861D+01 
0 . 4145443315 0D+01 
-0.758266 6042 6D+01 
0 . 31732538086D+01 


-0. 10248679150D+02 


0. 52914248276D+00 


THE  CC  MATRIX  (COLUMNS  1  -  4) 

0. 35961 167081D+01  0 . 15 1752 12 634D+01 

THE  CC  MATRIX  (COLUMNS  5  -  8) 

0.972 6292 6290D+01  0 . 20239432618D+01 
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